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Introduction 

In  this  work,  a  novel  imaging  technique  is  explored  that  uses  non-harmful  application  of  near  infrared  light  to  determine  the  properties 
of  tissue.  Using  this  technique,  known  as  near  infrared  (NIR)  tomography,  an  optical  fiber  placed  on  the  surface  of  the  region  of 
interest,  the  breast,  delivers  an  input  signal  while  other  optical  fibers  placed  at  different  locations  on  the  same  surface  detect  the  out- 
coming  photons,  which  have  propagated  through  the  volume  under  investigation.  The  intensity  and  path-length  distributions  of  the 
exiting  light  provide  information  about  the  optical  properties  of  the  transilluminated  tissue  using  a  model-based  interpretation  where 
photon  propagation  is  simulated  by  the  diffusion  theory.  Through  iterative  solution  to  match  the  theory  to  the  real  measurements, 
images  of  internal  absorption  and  scattering  coefficient  distribution  can  be  reconstructed.  Tissue  optical  properties  have  shown  to  be  a 
function  of  their  structure  and  more  importantly,  their  physiological  state.  It  is  these  qualities  that  provide  an  alternative  method  for 
the  detection  and  characterization  of  tumor  within  breast  tissue.  An  important  aspect  for  the  success  and  accuracy  of  this  method  is  the 
use  of  adequate  modeling  of  light  within  breast  tissue.  In  this  work,  we  have  used  a  Finite  Element  Model  of  the  Diffusion 
Approximation,  to  calculate  and  predict  light  transport.  Using  such  technique,  we  are  able  to  determine,  to  some  accuracy,  the 
distribution  of  internal  optical  properties  of  the  tissue  under  investigation.  A  major  limitation  in  the  development  of  near  infrared 
imaging  so  far  has  been  the  lack  of  accurate  and  fast  model-based  reconstruction  algorithm  for  this  modality. 

Body 

Since  the  award  of  this  project,  the  aims  have  not  changed.  In  brief  these  are:(l)  Validate,  improve  and  accelerate  a  three-dimensional 
finite  element  model  of  light  propagation  within  the  breast  tissue  as  well  as  the  image  reconstruction  algorithm  for  simultaneously 
calculating  the  internal  absorption  and  scattering  coefficients.  (2)  Investigate  the  benefits  and  limits  of  using  a-priori  information  from 
dual  modality  images,  for  example,  MRI  and  Near  Infrared  data.  (3)  Explore  the  benefits  and  limits  of  using  Fluorescent  contrast- 
agents,  which  will  target  specific  molecular  markers  of  the  tumor,  giving  rise  to  a-priori  information  regarding  tumor  location  as  well 
as  physiological  function  of  the  tissue.  (4)  Test  the  developed  algorithm  on  patient  data  and  compare  results  with  available 
mammograms  and  biopsy  results  to  evaluate  algorithm  accuracy 

Key  research  accomplishments 

In  the  initial  12  months  of  the  project,  much  has  been  achieved.  Namely,  (a)  Three-dimensional  reconstructed  images  of  patient  data 
have  been  achieved,  (b)  new  image  reconstruction  algorithms  have  been  developed  and  (c)  initial  work  has  begun  in  the  evaluation  of 
incorporating  a-priori  information  from  MRI  into  NIR  image  reconstruction. 


Reportable  outcomes 

Preliminary  Clinical  Results  of  Frequency  Domain  DOT 

Here,  we  present  sample  reconstructed  images  of  internal  absorption  and  scatter  from  2  sets  of  exams  where  data  was  collected  at 
multiple  wavelengths,  ranging  from  761  nm  to  826  nm  on  volunteers  recruited  at  Dartmouth-Hitchcock  Medical  Center.  The  NIR 
imaging  system  hardware,  developed  by  the 
research  team  at  Dartmouth,  is  shown  in  Figure  1. 

The  subject  lies  flat  on  the  measurement  bed 
Figure  1(a),  which  contains  a  single  opening  for 
the  breast.  The  breast  is  suspended  freely  through 
this  opening  and  the  optical  fibers,  Figure  1(b),  are 
brought  into  contact  with  the  breast  Figure  2.  The 
optical  fiber  arrangement  consists  of  a  total  of  48 
fibers,  16  fibers  in  3  separate  planes.  Within  each 
plane  the  fibers  are  arranged  equidistantly  (each 
separate  plane  is  10  mm  apart  from  one  another  in 
the  z  direction).  The  fibers  need  to  make  full 
contact  with  the  breast  for  high  quality  NIR 
measurements. 

For  the  NIR  exams,  an  attendant 
translates  the  fiber  optic  probes  into 
direct  contact  with  the  breast  at  the 
level  of  the  clinical  abnormality.  We 
obtained  three  tomographic  acquisitions 
centered  on  the  region  of  interest,  with 
contiguous  slices  above  and  below  the 
primary  plane  of  the  abnormality. 

Measurements  collected  from  a 
cylindrical  phantom  at  each  wavelength 
were  also  recorded  for  calibration 
purposes.  From  knowledge  of  the 
diameter  of  each  measurement  plane 
and  of  the  separation  between  planes, 


(a (b)l 
Figure  1.  (a)  The  exam  station  with  a  single  opening  for  the  breast,  (b)  Optical  fiber  probe 
arrangement  (16  equally  spaced  fibers  at  3  separate  levels,  each  level  is  separated  by  10  mm) 
currently  at  use.  The  fibers  can  be  opened  or  closed  to  suit  each  individual  subject  during  setup. 
The  breast  is  suspended  through  this  opening,  and  optical  fiber  probes  are  brought  into  contact 

with  the  tissue  surface. 


Figure  2.  Example  of  a  mesh  used  for  the  reconstruction  of  images  from  measured  clinical  data.  The 
information  about  data  collection  geometry  that  was  used  for  mesh  generation  included  slice  diameter, 
allowing  generation  of:  (a)  a  conical  shaped  mesh  used  for  light  field  calculation;  (b)  a  conical  shaped  mesh 

used  for  the  reconstruction  basis 
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we  constructed  a  conical  shaped  mesh,  Figure  2(a).  For  image  reconstruction,  each  data  set  was  calibrated  for  each  wavelength 
according  to  procedures  described 
elsewhere.  For  the  reconstruction 
basis  a  second  mesh  of  the  same 
geometry  was  used,  Figure  2(b). 

The  reconstruction  time  after  the 
initial  calibration  procedure  was 
approximately  10  min  per  iteration 
on  a  1.7-GHz  PC  with  2  Gigs  of 
RAM 

The  first  case,  (patient  A), 
presented  for  standard  screening 
mammography,  which  revealed  a 
subtle  nodular  density  and 
associated  architectural  distortion 
in  the  lateral  aspect  of  the  right 
breast.  Pathology  showed  an 
invasive  carcinoma  of  20  mm  size. 

3D  images  of  internal  absorption 
and  reduced  scattering  were 
recovered  simultaneously  from 
NIR  data  collected  at  each 
wavelength.  The  mesh  used  for  the 
calculation  of  the  Jacobian 
contained  8334  nodes 
corresponding  to  41623  linear 
tetrahedral  elements.  For 
reconstruction  basis,  a  second  mesh 
of  the  same  geometry  was  used,  but 
with  2521  nodes,  corresponding  to 
11575  linear  tetrahedral  elements. 

Images  were  reconstructed  at  four 
wavelengths  of  761,  785,  808  and  826  nm,  and  they  are  shown  in  Figure  3.  For  image  reconstruction  the  regularization  parameter  (k) 
used  was  initially  set  to  10,  and  was  allowed  to  decrease  by  a  factor  of  101/4  if  the  projection  error  had  decreased  with  respect  to 
previous  iteration.  Images  shown  are  those  at  the  10th  iteration.  Here  the  images  are  true  3D  reconstructions  and  coronal  slices  at  z  =  - 
60,  -45,  -30,  -15  and  0  mm  are  shown.  From  the  reconstructed  images  it  can  be  seen  that  an  anomaly  is  found  within  the  mid-plane  at 
approximately  9  o’clock  position.  The  anomaly  presents  as  an  absorption  variation  with  wavelength,  whereas  the  reduced  scatter  is 
almost  constant  over  all  reconstructed  wavelengths.  The  absorption  images  were  used  together  with  values  of  extinction  coefficients 
for  oxy  and  deoxy  hemoglobin  (Hb02  and  Hb  respectively)  for  the  calculation  of  Hb,  Hb02,  total  hemoglobin  (HbT)  and  oxygen 
saturation  (S02).  The  calculated  3D  maps  of  Hb,  Hb02,  total  HbT  and  S02  are  also  shown  in  Figure  3. 

From  the  calculated  values  of  Hb  it  is  seen 
that  the  anomaly  shows  a  peak  value  of  47.7 
pMole,  compared  to  a  background  of  32.87 
pM,  whereas  the  Hb02  image  shows  a  peak 
value  of  24.73  pM  at  a  location  on  the 
periphery  of  the  skin.  The  total  hemoglobin 
value  shows  a  peak  at  the  location  of  the 
anomaly,  with  a  value  of  65.28  pM.  S02 
values,  calculated  by  taking  the  ratio  of 
oxygenated  blood  and  total  blood  content, 
indicate  a  marked  decrease  at  the  location  of 
the  anomaly,  with  a  value  of  25.9%,  as 
compared  to  a  background  value  of  38.1 1%. 

These  data  are  summarized  in  Table  1 . 

Volunteer  B  presented  for  standard 
screening  mammography,  which  after  biopsy  revealed  a  benign  ductal  hyperplasia  region  within  the  right  breast.  From  the 
reconstructed  NIR  images  shown  in  Figure  4  an  anomaly  is  found  within  the  mid-plane,  off-center  near  the  nine  o’clock  position. 
Here,  the  anomaly  indicates  the  variation  of  absorption  with  wavelength,  with  a  peak  value  of  0.0055  mm”1  at  761  nm,  and  the  reduced 
scattering  images  show  a  peak  value  of  3.3  mm'1,  also  at  761  nm.  The  calculated  maps  of  Hb  shows  a  value  of  10.76  pMs  at  the 
location  of  the  anomaly,  compared  with  a  background  of  mean  value  of  7.7  pMs.  The  Hb02  image  shows  a  max  value  of  15.56  pMs, 


Table  1.  Values 
increase^ 

of  Hemoglobin  and  oxygen  saturation  maps  shown  in  Figures  15  and  16.  The  arrows  show  the 
b)  or  decrease  (4-)  of  values  within  the  region  of  interest  with  respect  to  the  background. 

Hb  Hb02  HbT  S02(%)  Tumor  Type 

(HM)  (uM)  (nM) 

Volunteer  A 
(anomaly) 

47.7  t 

24.73  + 

65.28  + 

25.9  4- 

Invasive  ductal  carcinoma 

Volunteer  A 
(background) 

32.87 

20.23 

53.11 

38.11 

Volunteer  B 
(anomaly) 

10.76  + 

15.56  ♦ 

22.56  + 

69.26  + 

Benign  ductal  hyperplasia 

Volunteer  B 
(background) 

7.71 

8.40 

16.11 

51.52 

Absorption 

761  nm  785  nm  808  nm  826  nm 


Reduced  Scatter 

761  nm  785  nm  808  nm  826  nm 


Hb  HbOz  HbT  S02 


Figure  3.  Reconstructed  images  of  absorption  (top  row)  and  reduced  scattering  (middle  row)  at  each  wavelength 
from  measured  volunteer  data,  (volunteer  A),  as  well  as  calculated  maps  of  blood  content;  each  slice  represents  a 
plane  through  the  mesh,  from  the  bottom  near  the  nipple  to  the  top  near  the  chest.  The  images  are  coronal  views  of 
the  cross  section  through  the  breast  at  the  tenth  iteration  at  the  waveleneths. 
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compared  with  a  mean  background 
value  of  8.4  pMs.  From  the  S02 
images,  the  anomaly  shows  a 
marked  increase  to  69.26%  in 
saturation,  compared  with  a  mean 
background  saturation  of  51.52% 

(Table  1). 

For  volunteer  A,  who  had  an 
invasive  ductal  carcinoma  within 
the  breast,  the  absorption  and 
scattering  increased  within  the 
region  of  interest.  The  calculated 
values  of  absorption  were  used 
together  with  extinction 
coefficients  of  Hb,  Hb02  to 
calculate  deoxy,  oxy  and  total 
hemoglobin  maps,  and  from  these 
data  oxygen  saturation  images 
were  generated.  While  the  total 
hemoglobin  level  within  the  region 
of  interest  shows  an  increase  with 
respect  to  background,  the  oxygen 
saturation  shows  a  marked 
decrease  within  the  same  region. 

This  trend  is  perhaps  expected  for 
malignant  tissue,  as  one  would 
expect  a  rise  in  blood  content,  due 
to  an  increase  in  blood  vessel 
density,  but  since  malignant  tumors 
are  more  active,  this  would  lead  to 
a  decrease  in  oxygen  saturation. 

For  volunteer  B,  who  had  a  benign  ductal  hyperplasia  within  the  breast,  the  absorption  and  scattering  also  increased  within  the  region 
of  interest.  For  a  benign  condition,  one  might  expect  a  rise  in  blood  content  due  to  an  increase  in  blood  vessel  density,  however  benign 
tumors  are  not  particularly  more  active,  suggesting  this  would  result  in  an  increase  in  oxygen  saturation  as  indicated. 

It  should  be  stressed  that  these  results  are  preliminary,  and  further  investigation  is  needed  before  one  can  claim  tumor  classification 
from  tissue  oxygenation  maps.  Further  work  is  in  progress  with  respect  to  building  a  larger  database  of  such  images,  thereby  allowing 
a  proper  and  definite  statistical  analysis  of  the  results. 

Spectral  imaging 

Near  Infra-Red  (NIR)  tomography  is  known  for  its  ability  to  distinguish  between  malignant  and  normal  tissue  based  on  high-contrast 
arising  from  intrinsic  processes  such  as  angiogenesis  and  hypoxia  in  tumors.  Absorption-based  parameters  can  be  recovered  such  as 
total  hemoglobin,  oxygen  saturation  and  water  fraction  as  well  as  scattering.  These  provide  measures  of  physiology  or 
pathophysiology  related  to  the  intra-vascular  and  extra-vascular  spaces  as  well  as  composition  of  the  tissue.  Typically  spectral 
tomography  approaches  must  make  use  of  only  sparse  data  at  discrete  wavelengths  instead  of  a  complete  spectrum,  resulting  in  an 
insufficient  measurement  set  for  a  unique  solution.  The  reconstruction  process  is  an  ill-posed  problem,  which  further  amplifies  this 
error  leading  to  image  artifacts  and  inaccuracy  in  quantifying  the  physiological  parameters  of  the  tissue.  The  current  method  illustrated 
here  reduces  the  noise  in  recovery  of  chromophores,  especially  significant  in  water  and  scattering,  by  eliminating  the  need  to 
reconstruct  the  absorption  and  scattering  coefficients  as  intermediate  endpoints.  Instead,  the  chromophore  concentrations  and  scatter 
parameters  are  reconstructed  directly  by  incorporating  the  known  linear  spectral  fit  and  scattering  relationship  with  wavelength  as 
constraints. 

The  method  uses  finite  elements  to  model  the  diffusion  equation  and  reconstructs  images  of  five  parameters:  oxyhemoglobin  (HbO), 
hemoglobin  (Hb),  water  (H20),  scatter  amplitude  (a)  and  scatter  power  (b),  with  no  assumptions  on  scatter  amplitude  (a)  or  scatter 
power  (b).  The  results  of  the  application  of  spectral  constraints,  in  addition  to  showing  a  reduction  in  noise  in  the  recovered 
chromophore  concentrations,  also  provide  evidence  of  reduced  sensitivity  to  noise  in  measurements  themselves. 

The  direct  spectral  reconstruction  method  has  been  tested  using  simulation  and  phantoms  studies.  As  an  illustration  here,  it  has  been 
applied  to  clinical  tomography  data  obtained  from  measurements  on  a  66  year  old  subject  with  a  5-10  mm  spiculated  mass  from 
mammography,  later  diagnosed  as  Invasive  Ductal  Carcinoma.  The  subject  underwent  the  NIR  exam  in  accordance  with  the 
Dartmouth  protocol,  and  written  consent  was  obtained  from  the  subject.  The  position  of  the  anomaly  was  provided  by  mammography 
and  found  to  be  at  11:30  o’clock  in  the  cranocaudal  view.  The  tumor  was  known  to  be  located  about  1cm  from  surface;  data  was 
collected  in  three  planes.  Reconstructed  images  are  shown  in  Figure  5,  for  the  plane  in  line  with  the  tumor.  The  localized  increase  in 
HbT  was  observable  and  the  contrast  available  was  1.7: 1.0  in  tumor  versus  background.  The  St02  image  showed  a  decrease  at  the 
location  of  the  tumor,  with  a  contrast  of  0.7:1.  The  corresponding  St02  image  with  the  conventional  technique  (not  shown  here)  was 
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Absorption 

761  nm  785  nm  808  nm  826  nm 


Reduced  Scatter 

761  nm  785  nm  808  nm  826  nm 


Figure  4  Reconstructed  images  of  absorption  (top  row)  and  reduced  scattering  (middle  row)  at  each  wavelength 
from  measured  volunteer  data,  (volunteer  B),  as  well  as  calculated  maps  of  blood  content;  each  slice  represents 
a  plane  through  the  mesh,  from  the  bottom  near  the  nipple  to  the  top  near  the  chest.  The  images  are  coronal 
views  of  the  cross  section  through  the  breast  at  the  tenth  iteration  at  the  wavelengths. 


noisier,  and  a  similar  decrease  was  not  observed.  The  water  image  showed 
heterogeneity  in  the  range  25  to  50%,  which  is  considerably  smaller  than  the  range 
of  0  to  100%  observed  with  the  old  method.  Tighter  data  ranges  are  similarly 
observed  in  the  scatter  images,  and  artifacts  in  the  separate  wavelength  technique 
have  been  completely  suppressed  in  the  spectral  approach. 

In  summary,  the  new  direct  spectral  reconstruction  technique  has  been  applied  to 
homogeneous  and  heterogeneous  data  and  the  results  have  consistently  shown 
reduced  artifacts  as  well  as  improved  quantitative  accuracy  in  recovered 
parameters.  With  this  new  technique,  higher  hemoglobin  content  at  site  of  tumor  is 
observed  in  the  clinical  case,  as  well  as  lower  oxygen  saturation.  The  water  and 
scattering  images  most  significantly  show  improvement  in  all  cases  by  reduction 
in  the  errors  due  to  decreased  artifacts. 

Incorporation  of  anatomical  data 

An  important  advance  has  been  the  culmination  of  several  years  of  effort  towards 
refining  an  MRI  compatible  NIR  tomography  system  which  works  comfortably 
with  a  breast  imaging  coil  to  sample  light  signals  transmitted  through  breast  tissue. 
The  system  is  shown  in  the  figure  6(a)  below,  with  the  ring  of  16  bifurcated  fibers 
appearing  in  the  positioning  array  Figure  2(b),  with  brass  loaded  springs,  which 
gently  push  the  fibers  against  the  breast.  A  series  of  tissue  simulating  phantoms 
were  created  with  gelatin,  and  a  three-layer  phantom  is  shown  in  the  figure  below, 
where  an  interior  layer  is  included  to  mimic  glandular-tissue  surrounded  by  adipose  tissue.  A  spherical  interior  object  is  also  included 
to  mimic  a  large  tumor.  These  images  have  been  used  in  a  series  of  basic  studies,  where  we  are  determining  the  optimal  way  to 
integrate  spatial  a  priori  information  into  the  NIR  tomography  reconstruction. 
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Figure  5.  Hbr  [mM],  St02  [%],  water  [%],  ‘a’  and  V 
images  are  shown  from  measurements  on  a  cancer 
subject  with  a  5-10  mm  Invasive  Ductal  Carcinoma  at 
the  1 1:30  o-clock  position. 


(a)  (b)  (c)  (d)  (e) 

Figure  6.  Image  of  the  patient  on  the  breast  coil  lying  on  the  MRI  bed  is  shown  in  (a).  An  image  of  the  breast  interface  is  shown  in  (b)  with  spring-loaded  optical  fibers 
arranged  in  a  circle.  A  gelatin  phantom  is  shown  in  (c)  with  varying  concentrations  of  scatterer,  absorber  and  copper  sulphate  solution,  allowing  simultaneous  MRI  and 
NIR  phantom  testing.  MRI  images  of  the  phantom  (d)  and  a  breast  (e)  are  shown  where  the  presence  of  the  array  can  be  seen  as  indents  in  the  phantom  image,  and  as 
exterior  fiducials  in  the  breast  image.  On  each  fiber,  a  sponge  ring  fiducial  is  placed,  allowing  accurate  imaging  of  the  fiber  tip  on  the  breast  tissue. 


Images  of  the  phantoms  used  in  this  system  are  shown  in  Figure  7,  with  segmentation  of  the  three  layers  in  the  phantom.  Without 
using  the  MRI  interior  structural  information,  the  image  reconstruction  of  this  phantom  would  appear  as  shown  in  Figure  7(b).  From 
the  phantom  data  we  note  that  it  is  possible  to  significantly  improve  the  fitting  of  optical  properties  using  a  priori  information.  This 
latter  approach  is  more  accurate  in  estimating  the  peak  value.  Interestingly,  if  the  regions  were  recovered  as  bulk  values,  reducing  the 
number  of  estimators  to  only  three,  this  approach  leads  to  incorrect  properties  for  the  middle  region.  Thus,  we  have  tentatively 
concluded  that  a  priori  information  is  better  applied  into  the  regularization  parameter,  than  through  segmentation  and  parameter 
reduction.  Detailed  analysis  of  the  tradeoffs  between  parameter  reduction  (regionizing  areas  to  be  represented  by  a  single  pair  of  //a 
and  nl  values)  versus  full  reconstruction  with  regionization  of  the  regularization  parameter  has  been  completed  on  phantom  data  sets. 
The  conclusions  to  date  indicate  that  regularization-based  constraints  are  much  more  stable  in  the  presence  of  noise,  and  that  the 
covariance  between  nodes  in  the  same  medium  is  a  key  factor  to  include  in  this  recovery.  Thus,  our  working  plan  is  to  utilize  interior 
information  as  carefully  chosen  variations  in  the  regularization  parameter. 


(a)  (b)  (c)  (d) 

Figure  7.  Cross  sectional  MRI  (a)  and  NIR  (b  and  c)  images  of  the  tissue  simulating  phantom  are  shown,  with  the  MRI  image  in  (a)  being  used  to  define  the  exterior 
and  interior  boundaries  between  regions.  In  (b)  our  standard  NIR  image  reconstruction  results  are  shown,  where  the  peak  absorption  coefficient  is  found  to  reach  only 
0.007  mm'1,  where  the  true  value  was  0.02  mm'1.  If  the  phantom  interior  regions  obtained  from  the  MR  image  are  used  to  vary  the  regularization  parameter,  the  peak 
value  increased  to  0.01 1  mm'1  (c),  however  if  covariance  between  nodes  is  also  incorporated  (d),  the  value  peaks  (0.01 8  mm'1)  nearer  the  true  level. 

In  Figure  8,  images  of  a  fatty /scattered  breast  are  shown  in  (a)  and  a  heterogeneously  dense  breast  in  (b),  where  in  the  latter 
breast,  it  is  possible  to  segment  the  two  major  tissue  types.  These  tissues  are  adipose  on  the  exterior  and  fibro-glandular  structure  on 
the  interior.  Without  using  these  segmented  regions,  the  typical  NIR  image  of  this  breast  would  be  as  shown  in  (c),  whereas  with  the 
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segmentation  of  the  regions,  the  two  tissue  types  can  be  fit  to  average  bulk  values,  which  agree  with  the  expected  numbers  for  these 
tissues  (see  figure  caption). 


(a)  (b)  (c)  (d) 

Figure  8.  Coronal  MRI  images  of  a  fatty/scattered  breast  are  shown  in  (a),  and  a  heterogeneously  dense  breast  in  (b)  in  the  breast  imaging  array,  with  the  fiber  fiducial 
markers  demarking  the  fiber  locations.  The  reconstructed  NIR  images  of  the  breast  shown  in  (b)  are  presented  in  (c)  without  a  priori  interior  information  and  (d)  with 
interior  spatial  information  applied.  The  values  of  the  absorption  coefficient  for  the  breast  tissue  types  in  the  fatty  and  glandular  regions  are  0.0035  mm'1  and  0.006 
mm’1. 

Conclusions 

The  project  is  continuing  to  reconstruct  NIR  images  of  clinical  data  to  build  the  database  of  3D  images.  Steps  have  also  been  taken  as 
part  of  this  project  to  achieve  other  aims  specified  by  this  project.  A  new  and  novel  image  reconstruction  algorithm  has  been 
developed  and  implemented  which  allows  the  direct  reconstruction  of  chromophore  concentration  and  scattering  functions  of  breast 
tissue.  This  not  only  allows  a  computationally  efficient  implementations,  but  also  since  all  the  data  from  all  wavelength  are  used 
coupled  together,  a  great  reduction  in  image  noise  and  artifact  is  seen.  Finally,  as  a  direct  result  of  the  2nd  aim,  the  project  is  beginning 
to  develop,  understand  and  implement  a  useful  algorithm  for  the  incorporation  of  a-priori  information,  not  only  in  phantom  studies, 
but  also  in  clinical  setups. 
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Abstract 

Near  infrared  tomography  (NIR)  is  a  novel  imaging  technique  that  can  be 
used  to  reconstruct  tissue  optical  properties  from  measurements  of  light 
propagation  through  tissue.  More  specifically  NIR  measurements  over  a  range 
of  wavelengths  can  be  used  to  obtain  internal  images  of  physiologic  parameters 
and  these  images  can  be  used  to  detect  and  characterize  breast  tumour.  To 
obtain  good  NIR  measurements,  it  is  essential  to  have  good  contact  between 
the  optical  fibres  and  the  breast  which  in-turn  results  in  the  deformation  of  the 
breast  due  to  the  soft  plasticity  of  the  tissue.  In  this  work,  a  tissue  deformation 
model  of  the  female  breast  is  presented  that  will  account  for  the  altered  shape 
of  the  breast  during  clinical  NIR  measurements.  Using  a  deformed  model  of 
a  breast,  simulated  NIR  data  were  generated  and  used  to  reconstruct  images 
of  tissue  absorption  and  reduced  scatter  using  several  assumptions  about  the 
imaging  domain.  Using  either  a  circular  or  irregular  2D  geometry  for  image 
reconstruction  produces  good  localization  of  the  absorbing  anomaly,  but  it 
leads  to  degradation  of  the  image  quality.  By  modifying  the  assumptions  about 
the  imaging  domain  to  a  3D  conical  model,  with  the  correct  diameter  at  the 
plane  of  NIR  measurement,  significantly  improves  the  quality  of  reconstructed 
images  and  helps  reduce  image  artefacts.  Finally,  assuming  a  non-deformed 
breast  shape  for  image  reconstruction  is  shown  to  lead  to  poor  quality  images 
since  the  geometry  of  the  breast  is  greatly  altered,  whereas  using  the  correct 
deformed  geometry  produces  the  best  images. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Near  infrared  (NIR)  tomography  is  an  emerging  alternative  imaging  method  used  to 
image  physiologic  parameters  of  biological  tissue  in  vivo  such  as  water  and  haemoglobin. 
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Measurements  of  light  propagation  (600-900  nm)  within  tissue  can  be  used  to  map  internal 
chromophore  concentrations  within  tissue.  Light  is  transmitted  through  tissue  using  multiple 
input  and  output  locations  on  the  surface  of  the  region  to  be  imaged,  similar  to  a  fan- 
beam  x-ray  computed  tomography  geometry,  but  using  optical  fibres  for  delivery  and  pick 
up  of  the  light  signals.  The  intensity  and  path-length  distributions  of  the  exiting  photons 
provide  information  about  the  optical  properties  of  the  transilluminated  tissue  using  a  model- 
based  interpretation  where  photon  propagation  is  simulated  by  diffusion  theory.  Using  these 
techniques  implemented  in  hardware  and  software,  NIR  optical  tomography  becomes  an 
inherently  three-dimensional  imaging  method  and  is  used  to  reconstruct  physiologically 
relevant  chromophore  distributions  from  the  region  under  investigation  (Eda  et  al  1999, 
Fantini  et  al  1999,  Boas  et  al  200 1 ,  Hebden  et  al  200 1 ,  Pogue  et  al  2001 ,  McBride  et  al  2002, 
Dehghani  et  al  2003a,  2003b,  2003c).  The  main  interest  in  this  study  lies  in  the  ability  to 
detect  and  characterize  tumours  within  the  female  breast  (Pogue  et  al  2001,  Brooksby  et  al 
2003).  Since  the  absorption  and  scattering  of  light  in  tissue  is  a  function  of  its  optical 
properties,  and  hence  its  physiological  state,  the  aim  is  to  obtain  images  of  internal  optical 
absorption  coefficient  //.a,  reduced  scattering  coefficient  /zs’  and  ultimately  images  of  total 
haemoglobin  and  oxygen  saturation  distributions.  These  images  should  in  principle  provide 
information  about  the  physiological  state  of  the  tissue  under  investigation  and  help  identify 
and  characterize  tumours  within  the  breast.  However,  in  recent  years  it  has  become  apparent 
that  the  external  shape  of  the  tissue  can  dramatically  impact  the  quality  of  the  reconstruction, 
and  this  is  largely  due  to  mismatch  between  the  model  prediction  and  the  actual  shape  of  the 
tissue  boundary.  In  this  study,  a  focused  effort  has  been  put  forward  to  examine  the  extent  to 
which  external  boundary  changes  will  degrade  image  quality  and,  more  importantly,  how  this 
could  be  corrected  with  appropriate  model-based  analysis. 

To  obtain  clinical  measurements  with  a  sufficient  signal-to-noise  ratio,  it  is  key  to  ensure 
that  good  contact  exists  between  the  optical  fibres  and  the  tissue.  More  specifically  to  our 
studies,  the  patient  lies  prone  on  the  measurement  bed,  which  contains  a  single  opening  for 
the  breast.  The  breast  is  suspended  freely  through  this  opening,  below  which  the  optical  fibres 
are  brought  into  contact  with  the  breast.  Normally,  the  optical  fibre  arrangement  consists  of 
a  total  of  48  fibres,  16  fibres  in  three  separate  planes.  However,  for  the  purpose  of  this  work 
only  one  plane  of  16  fibres  is  considered.  The  fibres  need  to  make  full  contact  with  the  breast 
for  a  good  and  adequate  NIR  measurement. 

The  breast  is  a  soft  tissue,  which  will  deform  and  alter  its  shape  on  the  application  of 
external  pressure.  The  amount  of  deformation  is  a  function  of  the  tissue  mechanical  properties 
and  the  amount  of  displacement  and  the  external  pressure  applied  by  the  optical  fibre  array.  In 
most  image  reconstruction  algorithms,  the  general  assumption  is  made  that  the  region  under 
investigation  is  a  uniform  circular  (2D)  or  conical  or  cylindrical  (3D)  domain.  Little  work  has 
been  done  to  evaluate  the  effect  of  any  incorrect  geometry  in  image  reconstruction.  For  simple 
symmetric  geometries,  for  example  a  3D  cylinder,  work  has  been  described  to  accurately 
calibrate  the  data  to  account  for  2D/3D  mismatch  (Hillman  et  al  2000).  More  recently,  work 
has  been  done  with  regard  to  neonatal  head  imaging  that  a  slight  change  to  the  geometry 
model,  both  in  terms  of  geometry  and  mesh,  will  result  into  large  change  in  modelled  data 
(Gibson  et  al  2003a). 

In  this  work,  the  effect  of  such  assumptions  is  investigated  by  creating  a  deformation  of 
the  breast  model  in  3D  to  generate  simulated  clinical  data.  Images  are  then  reconstructed  using 
various  assumptions  regarding  the  imaging  domain  including  a  circular  or  irregular  2D  models 
as  well  as  a  3D  conical  shaped  model  and  the  non-deformed  breast  model  to  evaluate  image 
quality.  An  important  part  of  this  modelling  effort  is  the  recent  advances  which  have  been  made 
in  modelling  soft  tissue  deformation  (Paulsen  et  al  1999,  Doyley  et  al  2000,  Van  Houten  et  al 
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2000).  Beyond  this,  it  is  also  becoming  established  that  soft  tissue  elastic  properties  can  be 
directly  measured  with  elastography  imaging  in  vivo ,  thereby  providing  the  key  information 
needed  for  predicting  the  deformation  of  tissue  under  force  or  displacement  conditions.  In  this 
study,  this  displacement  modelling  was  used  along  with  the  NIR  tomography  modelling 
to  create  a  comprehensive  three-dimensional  prediction  of  how  to  approach  appropriate 
modelling  of  deformed  tissue  that  is  being  imaged  with  NIR  tomography. 

2.  Theory 

2.1.  Breast  deformation  model 

Soft  tissues  exhibit  nonlinear  elastic  behaviour  (Fung  1993).  Nevertheless,  it  can  be  considered 
a  linear  elastic  material  in  situations  where  the  deforming  forces  produce  infinitesimal 
deformations  (i.e.  ^5%)  (Timoshenko  and  Goodier  1970).  For  the  purpose  of  this  study, 
the  breast  was  modelled  as  linear  isotropic  pseudo-incompressible  (i.e.  Poisson’s  ratio  (v)  = 
0.495  (Fung  1993)).  Under  these  assumptions  and  ignoring  internal  body  force,  the  governing 
elasticity  equations  for  quasi-static  deformation  are  given  by  Timoshenko  and  Goodier  (1970) 
and  Fung  (1993) 

(X  +  pi)V(V -u)  +  piV2u  =  0  (1) 

for  internal  nodes  in  domain  ft,  and 

((k  +  /i)V(V  •  u)  +  pxV2u)  •  h  =  h  (2) 

for  nodes  on  the  boundary  5ft. 

Here  h  represents  a  unit  vector  directed  outwards  from  ft,  and  h  represents  the  traction  on 
the  surface  or  boundary  of  the  breast.  Note  that  u  represents  the  displacement  components  in 
all  coordinate  directions,  and  pi  and  y  are  Lame’s  elastic  constants.  For  an  isotropic  medium 
these  constants  are  related  to  the  more  familiar  Young’s  modulus  ( E)  and  the  Poisson’s  ratio 

00  by 

E  vE  (3) 

M  2(1  +v)  (1  +  v)(l  —  2v) 

The  first  Lame’s  elastic  constant  (i.e.  pi)  is  generally  known  as  shear  modulus. 

It  should  be  stated  that  in  this  study,  the  breast  was  assumed  to  be  traction  free  (i.e.  no 
other  forces  are  associated)  and  internal  body  forces  were  neglected,  thus  the  problem  is  solved 
by  imposing  a  prescribed  displacement  (i.e.  induced  when  the  optical  fibres  are  coupled  to  the 
breast)  as  described  in  Doyley  et  al  (1999). 

2.2.  Light  propagation  model 

Under  the  assumption  that  scattering  dominates  absorption  for  NIR  light  in  tissue,  the 
Boltzmann  transport  equation  can  be  simplified  to  the  diffusion  approximation,  which  in 
the  frequency  domain  is  given  by 

-V  •  DVO(r,  co)  +  (p,a  +  <t>(r,  co)  =  qo(n  o>)  (4) 

where  q0(r ,  co)  is  an  isotropic  source,  3>(r,  co)  is  the  photon  fluence  rate  at  position  r  and 
D  =  -  is  the  diffusion  coefficient.  We  use  the  Robin  (Type  III)  boundary  condition: 

iyHa+r-s) 

4>(y)  +  —h  -  V3>(y)  =  0 
a 


(5) 
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where  a  is  a  term  which  incorporates  reflection  as  a  result  of  refractive  index  mismatch 
(Schweiger  et  al  1995,  Dehghani  et  al  2003c)  at  the  boundary,  and  n  is  the  outward  pointing 
normal  to  the  boundary  (8Q)  at  y. 

We  assume  that  the  data  are  represented  by  a  nonlinear  operator  y*  =  D],  where 

our  data  y *  are  a  complex  vector  having  a  real  and  imaginary  components,  which  are  mapped 
to  log  amplitude  and  phase  shift  in  measurement.  Then  the  image  reconstruction  method  is 
posed  as  a  solution  to  the  following  expression: 

(Ao,  D)  =  argmin  ||(;y*  -  F(jia,  D))||  (6) 

where  ||  •  ||  is  the  weighted  L2-norm,  representing  the  square  root  of  the  sum  of  the  squared 
elements,  fta  is  a  vector  of  the  absorption  coefficients  and  D  is  a  vector  of  the  diffusion 
coefficients.  The  magnitude  of  this  is  sometimes  referred  to  as  the  projection  error  and 
provides  a  value  for  determining  the  convergence  of  the  iterative  reconstruction  algorithm. 

In  this  study,  a  finite-element  method  (FEM)  is  used  as  a  general  and  flexible  method  for 
solving  the  forward  problem  in  arbitrary  geometries  (Arridge  et  al  1993,  Jiang  et  al  1996).  In 
the  inverse  problem,  where  the  goal  is  to  recover  internal  optical  property  distributions  from 
boundary  measurements,  it  is  assumed  that  fia( r)  and  D( r)  are  expressed  in  a  basis  with  a 
limited  number  of  dimensions  (less  than  the  dimension  of  the  finite  element  system  matrices). 
A  number  of  different  strategies  for  defining  reconstruction  bases  are  possible;  in  this  paper  a 
linear  pixel  basis  (Schweiger  and  Arridge  1999)  is  used.  To  find  (£fl,  D)  in  equation  (6)  we 
have  used  a  Levenberg-Marquardt  algorithm,  where  we  repeatedly  solve: 

a  =  JT  (J  JT  +  pl)~xb  (7) 

where  b  is  the  data  vector,  b  =  [y*  -  F[p,a ,  D]]r,  a  is  the  solution  update  vector, 
a  =  [ 8D\8jxa ].  Here,  p  is  the  regularization  factor  and  J  is  the  Jacobian  matrix  for  our 
model  which  is  calculated  using  the  Adjoint  method  (Arridge  and  Schweiger  1995,  Dehghani 
et  al  2003b,  2003c). 

3.  Methods 

3.1.  Breast  deformation 

A  volume  mesh  of  a  female  breast  of  a  volunteer  was  created  from  surface  image  data  that  was 
acquired  using  a  3D  surface  camera  (Rainbow  3D  Camera,  Genex  Technologies,  Kensington 
MD).  The  3D  camera  projects  structural  illumination  patterns  onto  the  object  and  calculates 
3D  surface  profile  described  by  over  300000  data  points  (Geng  1996,  Galdino  et  al  2002).  A 
volume  mesh  is  then  generated  using  the  Delaunay  algorithm.  The  mesh  has  a  geometry  of 
130  x  136  x  60  mm,  figure  1(a),  and  contained  15  501  nodes  corresponding  to  61  171  linear 
tetrahedral  elements.  The  diameter  of  the  breast  at  its  mid-plane  where  optical  fibre  array  will 
be  applied  is  approximately  88  mm.  To  calculate  the  deformation  due  to  16  equally  spaced 
optical  fibres  being  applied  at  the  mid-plane  of  the  breast,  i.e.  z  =  —30  mm,  it  is  assumed  that 
each  optical  fibre  pushed  the  breast  so  that  the  final  breast  diameter  at  z  =  —30  mm  is  70  mm, 
and  that  the  diameter  of  each  optical  fibre  is  6  mm.  The  modelled  elastic  properties  of  tissue, 
equation  (3),  were  assumed  as  isotropic  and  homogenous  with  Young’s  modulus  of  20  kPa 
(Krouskop  et  al  1998)  and  Poisson’s  ratio  of  0.495  (Fung  1993).  Further,  it  was  assumed  that 
the  topmost  part  of  the  mesh,  i.e.  z  =  0  mm  was  not  allowed  to  move  since  it  is  connected 
to  the  chest.  Using  this  applied  displacement  as  a  boundary  condition,  the  displacement  at 
all  nodes  due  to  the  application  of  the  optical  fibres  was  calculated  and  a  deformed  mesh  was 
created,  figure  1(b). 
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Figure  1.  Volume  mesh  of  the  (a)  normal  suspended  breast  and  (b)  the  deformed  mesh  after  the 
application  of  the  optical  fibre  array. 


32.  Simulation  of  data  from  deformed  breast 

In  order  to  accurately  simulate  the  clinical  settings,  two  localized  anomaly  regions  were  placed 
within  the  mid-plane  of  the  normal  breast,  both  at  z  =  —30  mm,  figure  2(a).  First  was  an 
absorbing  anomaly,  19  mm  from  the  surface  and  a  radius  10  mm  with  a  coefficient  value  of 
0.03  mm-1.  Second  was  a  reduced  scattering  anomaly,  also  19  mm  from  the  surface  but  with 
a  radius  of  7.5  mm  and  a  value  of  3  mm-1.  These  anomalies  were  chosen  as  they  represent 
the  size  and  contrast  we  aim  to  image  and  so  that  absorption  and  scatter  separation  can  also  be 
investigated.  The  background  optical  properties  were  modelled  with  an  absorption  coefficient 
of  0.01  mm"1  and  a  reduced  scatter  coefficient  of  1.0  mm-1.  For  the  deformed  breast  model 
the  same  optical  properties  of  the  normal  breast  were  assumed  (the  anomalies  were  also 
assumed  to  have  the  same  elastic  properties  of  the  background).  Once  the  deformation  of  the 
model  was  calculated,  it  was  assumed  that  the  anomalies  were  free  to  move,  depending  on 
the  elastic  properties  of  the  breast,  figure  2(b).  It  is  interesting  to  note  that  the  displacement  of 
the  anomalies  is  not  so  much  dependent  on  the  mesh  density  of  the  breast,  but  more  dependent 
on  the  mechanical  properties,  applied  pressure  and  non-symmetric  nature  of  the  breast. 

Using  the  deformed  mesh,  together  with  the  displaced  anomalies,  data  were  simulated  for 
16  equally  circularly  spaced  optical  fibres  placed  at  z  =  —30  mm.  Amplitude  and  phase  data 
were  simulated  at  100  MHz,  and  1%  noise  was  added.  These  data  were  then  used  as  simulated 
patient  data  in  the  following  sections.  In  addition,  to  allow  data  calibration  as  done  using 
clinical  data,  and  discussed  elsewhere  (Dehghani  et  al  2003c,  McBride  et  al  2003),  the  data 
were  simulated  for  16  equally  spaced  optical  fibres  placed  in  a  circle  around  the  mid-plane  of  a 
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z  =  -50  mm  z  =  -40  mm  z  =  -30mm  z  =  -20  mm  z  =  -10  mm 


Absorption 


z  =  -50  mm  z  -  -40  mm  z  =  -30  mm  z  =  -20  mm  z  =  -10  mm 


Reduced  Scatter 

(b) 


Figure  2.  2D  coronal  slices  through  (a)  the  normal  breast  mesh  and  (b)  the  deformed  breast  mesh, 
showing  the  position  of  the  anomalies.  The  most  right-hand  slice  is  near  the  chest  while  the  most 
left-hand  slice  is  near  the  nipple. 
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Figure  3.  2D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast 
model  simulated  data.  The  mesh  used  for  the  reconstruction  is  a  circular  mesh  whose  diameter  is 
the  same  as  the  optical  fibre  array  diameter  used  to  deform  and  simulate  data  NIR  from  the  breast 
model. 


non-deformable  cylindrical  model  with  a  known  homogenous  background  optical  properties. 
The  same  calibration  data  file  has  been  used  for  all  presented  reconstructions. 

3.3.  Image  reconstruction  using  2D  meshes 

Using  the  computed  data  for  the  deformed  mesh  and  the  reference  phantom,  as  simulated 
measurements  with  1  %  random  noise  added,  the  data  were  calibrated  (Dehghani  et  al  2003c, 
McBride  et  al  2003)  and  images  were  reconstructed  using  two  different  2D  meshes.  Assuming 
that  the  imaging  domain  was  a  circular  domain  with  a  radius  equal  to  the  radius  of  the  optical 
fibre  array,  i.e.  35  mm,  images  were  reconstructed  and  the  results  are  shown  in  figure  3.  The 
circular  mesh  used  contained  1785  nodes  corresponding  to  3418  linear  triangular  elements. 
For  the  reconstruction  basis  (basis  used  for  the  update  of  the  optical  parameters),  a  20  x  20 
regular  grid  (piecewise  linear)  was  used  and  the  initial  regularization  parameter  was  set  to  100. 
Images  shown  here,  and  all  following  sections  are  at  iteration  level  where  the  projection  error 
was  within  5%  of  the  previous  iteration.  In  the  2D  reconstruction  cases,  this  projection  error 
was  reached  typically  in  the  7th  iteration. 

It  is  evident  from  figure  1(b)  that  the  breast  at  the  plane  of  measurement,  i.e.  z  =  —30  mm, 
was  no  longer  circular.  In  order  to  evaluate  if  a  more  correct  a  priori  information  regarding 
the  2D  boundary  at  the  plane  of  measurement  would  be  useful,  a  2D  irregular  mesh,  with  a 
boundary  whose  profile  is  the  same  as  the  boundary  at  z  =  —30  mm  for  the  deformed  mesh, 
was  created  and  used  for  image  reconstruction.  The  resulting  images  reconstructed  on  this 
mesh  are  shown  in  figure  4.  The  irregular  mesh  used  contained  2405  nodes  corresponding 
to  4619  linear  triangular  elements.  For  the  reconstruction  basis,  a  20  x  20  regular  grid  was 
used  and  the  initial  regularization  parameter  was  set  to  100.  Thus  in  this  case,  the  irregular 
boundary  was  matched,  but  still  the  image  reconstruction  was  approximated  by  a  2D  forward 
solution. 

3.4.  Image  reconstruction  using  3D  conical  mesh 

In  the  next  section,  we  examined  the  improvement  if  a  true  3D  forward  calculation  was 
assumed,  but  with  a  regular  geometry.  In  a  previous  work,  it  was  hypothesized  that  given 
a  set  of  3D  patient  data,  those  images  could  be  reconstructed  with  reasonable  accuracy  if 
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Absorption  Reduced  Scatter 


Figure  4.  2D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast 
model  simulated  data.  The  mesh  used  for  the  reconstruction  is  an  irregular  whose  boundary  is  the 
same  as  the  boundary  of  the  deformed  breast  mesh  at  z  =  30  mm. 


Figure  5.  The  conical  shaped  volume  mesh.  The  geometry  of  the  mesh  is  such  that  the  diameter  at 
the  mid-plane  corresponds  to  the  diameter  of  the  optical  fibre  array  used  to  deform  and  simulate  NIR 
data  of  the  breast  model.  The  diameter  at  20  mm  above  and  20  mm  below  the  NIR  measurement 
plane  of  the  deformed  breast  were  also  used  to  set  the  upper  and  lower  bound  of  the  cone  model. 


the  breast  was  assumed  to  be  a  conical  shaped  model  (Dehghani  et  al  2003c).  To  create  a 
conical  mesh,  the  information  about  the  diameter  of  the  measurement  fibre  array  as  well  as 
the  diameter  of  the  breast  at  20  mm  above  and  20  mm  below  the  measurement  plane  was 
used,  which  would  include  all  boundary  information  other  than  effects  due  to  tissue  bulging. 
Following  this  procedure,  a  conical  shaped  volume  mesh  was  created,  as  shown  in  figure  5, 
which  shows  the  upper  chest  wall  and  the  confined  partial-conical  breast  shape  as  a  regular 
geometry  approximation  to  the  circularly  compressed  breast  shape.  The  mesh  contained  9332 
nodes  corresponding  to  42281  linear  tetrahedral  elements  and  was  created  using  NETGEN 
(Schoberl).  Using  this  mesh  and  the  simulated  data  for  the  deformed  breast  mesh,  images 
were  reconstructed  and  the  results  are  shown  in  figure  6. 

For  the  reconstruction  basis,  a  20  x  20  x  10  regular  grid  was  used  and  the  initial 
regularization  parameter  was  set  to  100.  The  images  shown  here  and  all  through  the  following 
section  are  at  iteration  level  where  the  projection  error  was  within  5%  of  the  previous  iteration. 
These  presented  from  the  3D  cases  are  from  the  14th  iteration. 


3.5.  Image  reconstruction  using  the  normal  and  deformed  breast  mesh 

Assuming  that  correct  information  regarding  the  breast  is  available,  before  the  application 
of  the  optical  fibre  array,  the  normal,  non-deformed  mesh  was  used  as  shown  in  figure  1(a), 
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Figure  6.  3D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast 
model  simulated  data.  The  mesh  used  for  the  reconstruction  is  the  mesh  shown  in  figure  5.  The 
most  right-hand  slice  is  near  the  chest  while  the  most  left-hand  slice  is  near  the  nipple. 


as  the  reconstruction  mesh  for  the  images  shown  in  figure  7.  The  reconstruction  basis  and 
parameters  used  were  the  same  as  for  the  conical  model  already  described.  Finally,  to  show 
the  best  possible  reconstruction,  the  deformed  mesh,  figure  1(b)  was  used  for  reconstruction 
of  images  and  these  are  shown  in  figure  8. 


4.  Results 

The  deformed  breast  mesh,  as  a  result  of  applying  a  circular  optical  fibre  array  at  z  =  30  mm, 
is  shown  in  figure  1(b).  It  is  evident  that  the  mesh  has  been  compressed  at  the  plane  of  applied 
displacement  with  slight  bulging  in-between  each  optical  fibre.  The  breast  mesh  has  also 
extended  in  the  z  direction  near  the  nipple  to  maintain  a  constant  volume  as  the  nodes  on  the 
chest  wall  are  assumed  fixed.  Although  the  mechanical  properties  of  the  tissue  are  assumed 
constant,  it  is  interesting  to  note  that  the  modelled  anomalies  have  also  been  displaced  due  to 
the  deformation,  figures  2(a)  and  (b).  Both  the  anomalies  have  been  reduced  in  diameter,  but 
have  extended  in  the  z  direction  as  a  result  of  the  applied  deformation. 

Two-dimensional  images  reconstructed  using  the  simulated  deformed  breast  data  are 
shown  in  figures  3  and  4.  In  both  cases  where  a  circular  or  an  irregular  boundary  was 
used,  good  images  were  recovered  in  terms  of  localization  of  the  absorption  anomaly  (within 
1.2  mm  for  the  circular  model,  and  2.0  mm  for  the  irregular  model).  The  recovery  of  the 
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Figure  7.  3D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast 
model  simulated  data.  The  mesh  used  for  the  reconstruction  is  the  mesh  shown  in  figure  1(a).  The 
most  right-hand  slice  is  near  the  chest  while  the  most  left-hand  slice  is  near  the  nipple. 


reduced  scattering  anomaly  was  not  as  good  as  the  absorption,  and  both  sets  of  images  contain 
large  boundary  artefacts.  The  background  values  of  both  the  absorption  and  reduced  scatter 
were  close  to  the  expected  value,  however,  the  target  absorption  value  is  less  than  50%  than 
expected. 

Three-dimensional  images  reconstructed  assuming  a  conical  shaped  breast  are  shown  in 
figure  6.  In  that  case,  both  the  absorption  and  the  scattering  anomalies  were  recovered  with 
good  localization  (within  2.00  mm  for  both  the  absorption  and  scattering  objects)  and  with 
little  boundary  artefacts.  Since  the  measured  data  were  about  the  mid-plane  of  the  breast, 
and  the  modelled  cone  geometry  also  had  the  optical  fibres  around  its  mid-plane,  the  objects 
recovered  are  centred  at  z  =  -30  mm.  Here  the  background  values  of  absorption  and  scatter 
were  a  little  lower  than  expected  values,  however,  the  absorption  anomaly  had  a  maximum 
value  of  0.014  mm-1  and  the  scattering  object  had  a  maximum  value  of  2  mm-1. 

Images  reconstructed  using  the  assumption  of  correct  non-deformed  breast  geometry  are 
shown  in  figure  7.  In  this  case,  although  the  absorbing  anomaly  has  been  reconstructed,  the 
reduced  scattering  images  contain  large  artefacts.  The  background  absorption  and  scatter 
values  are  about  0.006  mm-1  and  0.5  mm-1,  respectively,  while  the  maximum  absorption  for 
the  anomaly  is  0.014  mm’1.  Finally,  images  reconstructed  assuming  correct  3D  deformed 
boundary  are  shown  in  figure  8.  Here,  both  the  absorbing  and  scattering  regions  have  been 
recovered  within  less  than  1 .00  mm  of  expected  location,  giving  great  localization  accuracy. 
There  exists  some  cross  talk  between  the  absorbing  and  scattering  anomalies.  The  background 
absorption  and  scatter  values  are  about  0.01  mm"1  and  0.98  mm"1,  respectively,  while  the 
maximum  absorption  and  scatter  values  for  the  anomaly  is  0.013  mm"1  and  1.3  mm"1, 
respectively. 
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Figure  8.  3D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast 
model  simulated  data.  The  mesh  used  for  the  reconstruction  is  the  mesh  shown  in  figure  1(b).  The 
most  right-hand  slice  is  near  the  chest  while  the  most  left-hand  slice  is  near  the  nipple. 


5.  Discussions 

The  deformation  of  the  breast  due  to  the  application  of  16  circular  equally  spaced  optical 
fibres  has  been  modelled  using  a  computational  mechanic  model.  Using  this  deformed 
model  which  also  contains  a  single  absorbing  and  a  single  scattering  perturbations,  data  were 
simulated  and  1  %  noise  was  added  to  more  adequately  represent  clinical  measurements.  Using 
these  simulated  data,  reconstructed  images  of  absorption  and  reduced  scatter  were  generated 
simultaneously  using  various  different  geometries  of  the  mesh  upon  which  the  reconstruction 
would  take  place.  In  the  first  case,  it  was  assumed  that  the  model  used  for  reconstruction  was 
circular  with  a  radius  equal  to  the  radius  of  the  optical  fibre  array  after  breast  compression, 
which  was  35  mm.  The  reconstructed  images  have  shown  good  accuracy  in  the  localization  of 
the  absorbing  object  with  a  maximum  absorption  value  of  0.014  mm-1.  The  recovery  of  the 
scattering  objects  is  much  poorer  than  the  absorbing  anomaly,  which  can  be  attributed,  perhaps, 
to  its  small  size  (15  mm  diameter).  There  exist  boundary  artefacts  in  both  the  absorbing  and 
scattering  images.  It  is  possible  to  improve  the  quality  of  the  reconstructed  2D  images  by 
using  several  methods,  including  the  use  of  various  regularization  parameters  and/or  use  of 
a  priori  information.  A  good  method  to  improve  the  reconstructed  images  in  this  case  would 
be  to  segment  the  reconstructed  images  shown  in  figure  3,  and  use  these  as  a  priori  information 
to  more  accurately  obtain  quantitative  results.  The  application  of  such  algorithm  is  discussed 
elsewhere,  and  has  been  found  to  provide  superior  results  (Srinivasan  et  al  2004).  In  the  case 
of  knowing  and  using  correct  2D  boundary  information  in  image  reconstruction,  figure  4,  no 
useful  improvement  was  seen.  Although  correct  information  was  added  to  the  reconstruction 
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with  regard  to  the  boundary,  the  image  reconstruction  algorithm  was  still  constrained  to 
be  2D.  While  it  has  been  shown  many  times  that  the  2D  geometry  can  be  used  to  recover 
accurate  image  data  from  NIR  tomography  with  a  circular  boundary  (Pogue  and  McBride 
1999),  it  is  not  surprising,  perhaps,  to  find  that  this  simplification  does  not  continue  to  work 
as  the  exterior  boundary  becomes  more  complex  in  shape.  The  reconstruction  mesh  is  no 
longer  symmetric  and  contains  many  irregular  edges,  which  gives  rise  to  high  frequency 
noise  at  the  boundary,  and  those  changes  occur  symmetrically  around  each  optical  fibre 
making  the  calculation  even  more  dependent  upon  the  3D  geometry  out  of  the  plane  of 
imaging. 

Next  for  the  case  where  the  imaging  domain  was  assumed  to  be  conical  in  shape,  the 
reconstructed  images  were  shown  in  figure  6.  Here  it  was  assumed  that  the  mid-plane  of 
the  cone  had  a  radius  equal  to  the  radius  of  the  optical  fibre  array  after  breast  compression. 
Further,  it  was  assumed  that  the  measurements  of  the  diameter  of  the  breast  20  mm  above 
and  20  mm  below  the  measurement  plane  were  known.  Using  these  data,  a  conical  shaped 
mesh  could  be  created  for  any  patient  exam  as  shown  in  figure  5.  The  reconstructed  images 
had  recovered  both  the  absorbing  and  scattering  anomalies  with  reasonable  success.  The 
calculated  background  values  for  absorption  and  scatter  were  0.007  mm"*1  and  0.94  miti"1, 
respectively.  There  exists  a  small  cross  talk  between  the  absorption  and  reduced  scatter  images, 
plus  artefacts  at  the  boundary  of  the  images.  The  qualitative  and  quantitative  accuracy  of  the 
conical  shaped  model  reconstruction  are  better  than  the  2D  circular  case,  simply  because  a 
more  accurate  3D  model,  rather  than  a  2D  model  is  used. 

In  the  case  where  it  is  assumed  that  the  exact  breast  geometry  was  known,  but  have 
ignored  information  regarding  the  compression  due  to  the  optical  fibres,  the  resulting  images 
are  shown  in  figure  7.  Here,  although  the  absorption  anomaly  was  recovered  with  relatively 
good  accuracy,  the  reduced  scattering  image  contains  a  large  artefact.  Furthermore,  the 
calculated  background  values  for  absorption  and  scatter  are  0.006  mm"1  and  0.5  mm”1, 
respectively,  which  are  much  lower  than  expected.  These  images  are,  perhaps,  not  as  accurate 
and  useful  as  the  conical  geometry  case,  since,  although  there  is  a  more  accurate  3D  model, 
the  diameter  of  the  breast  was  not  correct,  particularly  within  the  measurement  plane,  i.e.  z  — 
—30  mm.  Finally,  for  the  case  of  using  a  geometrically  correct  deformed  mesh  to  reconstruct, 
the  images  are  shown  in  figure  8.  Here,  the  reconstructed  images  of  absorption  and  scatter 
were  found  with  superior  localization  accuracy.  The  scattering  object  has  also  been  recovered 
in  this  case.  However,  the  quantitative  accuracy  of  the  recovered  anomalies  is  not  as  good 
with  maximum  values  for  absorption  and  scattering  objects  of  0.013  mm”1  and  1.3  mm”1, 
respectively.  This  low  quantitative  accuracy  has  been  reported  elsewhere  and  is  a  common 
problem  in  3D  NIR  imaging  algorithms  (Dehghani  et  al  2003b,  Gibson  et  al  2003b),  which 
may  be  solved  using  of  multistage  reconstruction  algorithms  (Srinivasan  et  al  2004)  or  with 
inclusion  of  a  priori  data  (Brooksby  et  al  2003).  Finally,  using  the  correct  model  with  correct 
information  regarding  the  deformed  boundary  has  produced  images  with  very  little  or  no 
artefacts. 


6.  Conclusions 

In  this  work,  we  have  presented  a  breast  deformation  model  that  would  account  for  the 
change  in  shape  and  geometry  of  the  breast  due  to  the  application  of  a  circular  array  of 
optical  fibres  for  NIR  measurements.  The  proposed  model,  at  present,  assumes  that  the  breast 
has  homogenous  mechanical  properties,  which,  although  not  accurate,  provide  a  good  initial 
estimate  for  modelling  of  any  deformation.  Future  studies  may  focus  on  using  spatially 
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distributed  mechanical  properties,  as  these  can  be  imaged  with  good  spatial  resolution  with 
MR  elastography  (Paulsen  et  al  1999,  Van  Houten  et  al  2000,  Weaver  et  al  2001,  Doyley  et  al 
2003).  Also,  in  this  work,  we  have  deformed  the  breast  more  than  adequately  needed  by 
19  mm,  whereas  in  a  real  clinical  setting,  the  deformation  will  be  of  lower  magnitude,  typically 
about  5-10  mm.  In  practice  for  our  imaging  exams,  we  currently  attempt  to  deform  the  breast 
as  minimally  as  possible,  however,  there  are  strategic  benefits  to  compressing  in  that  the  signal 
transmitted  can  be  higher,  and  there  can  be  pressure-induced  changes  which  might  provide 
meaningful  contrast  about  what  the  tissue  is  composed  of  (Jiang  et  al  2003).  Nevertheless  we 
have  chosen  such  a  large  deformation  to  provide  us  with  a  worst-case  scenario  for  NIR  image 
reconstruction,  and  to  examine  if  it  is  feasible  to  exploit  larger  magnitude  deformations 
in  clinical  studies  while  not  compromising  the  integrity  of  the  data  collected  from  the 
breast. 

Using  simulated  data  from  the  deformed  breast  that  also  contains  optical  anomalies,  we 
have  reconstructed  images  assuming  various  reconstruction  geometries.  In  the  first  case 
of  assuming  that  the  reconstruction  mesh  is  either  a  circular  or  irregular  2D  mesh,  the 
reconstructed  images  have  shown  good  accuracy  in  recovering  the  location  of  the  absorbing 
anomaly.  The  2D  reconstructed  scattering  images  are  not  as  good  as  expected,  which  is  perhaps 
due  to  the  fact  that  the  scattering  object  had  a  small  diameter  of  15  mm  and  a  contrast  which 
was  three  times  the  background,  and  thus  the  image  being  dominated  by  boundary  artefacts. 
Also  it  is  evident  from  the  reconstructed  images  that  there  exist  large  boundary  artefacts  which 
are  due  to  the  incorrect  model.  However,  both  the  quality  and  quantitative  accuracy  of  2D 
image  reconstruction  will  be  improved  by  the  incorporation  of  a  priori  information  as  well  as 
a  multi-step  image  reconstruction  where  regions  of  interest  can  be  identified  and  isolated  for 
regional  reconstruction. 

In  the  case  of  3D  image  reconstruction,  assumptions  regarding  the  breast  not  being 
deformed  might  be  used  and  the  reconstructed  images  provide  relatively  good  absorption 
distributions,  but  the  images  contain  large  artefacts  which  will  likely  confound  their 
interpretation.  The  reason  for  this  can  be  explained  by  considering  the  overall  shape  of 
the  breast  before  and  after  deformation.  The  non-deformed  mesh  at  the  plane  of  measurement 
has  radius  that  is  19  mm  larger  than  the  deformed  mesh.  Additionally,  as  evident  from 
figures  1(a)  and  1(b),  it  is  found  that  the  breast  expands  above,  below  and  between  each  optical 
fibre  once  compressed.  This  large  change  in  shape  of  the  breast  contributes  to  the  large  artefacts 
seen  within  the  reconstruction.  However,  if  one  assumes  that  the  diameter  of  the  breast  is 
known  within  the  measurement  plane,  as  well  as  an  approximate  diameter  above  and  below 
it,  it  is  possible  to  create  a  pseudo-3D  conical  shaped  mesh  for  image  reconstruction.  Images 
reconstructed  using  this  3D  conical  mesh  have  shown  a  better  accuracy  in  recovering  both 
the  absorbing  and  scattering  anomalies.  Finally,  if  accurate  knowledge  regarding  the  breast 
deformation  is  available,  images  of  optical  properties  can  be  reconstructed  which  localize  both 
anomalies  with  great  accuracy,  and  also  contain  little  or  no  artefact  which  otherwise  would 
arise  from  model  mismatch. 

The  goal  of  this  work  is  to  incorporate  this  new  model  of  breast  deformation  with  more 
accurate  information  regarding  the  mechanical  properties  of  the  breast  to  improve  the  NIR 
image  reconstruction.  The  mechanical  property  information  is  readily  available  from  other 
imaging  modalities,  and  the  synthesis  of  this  information  may  provide  fundamentally  new 
information  about  breast  physiologic  response  to  pressure,  and/or  breast  pathology  response 
to  pressure.  An  accurate  model  of  breast  deformation  should  in  principle  allow  us  to  create 
patient  specific  models  and  meshes,  which  would  in-tum  provide  more  clinically  useful  data. 
Work  is  in  progress  to  validate  and  further  improve  the  deformation  model  with  application  to 
the  female  breast  imaging. 
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Abstract:  In  NIR  tomography  of  the  breast,  good  contact  is  needed  by  the  fibers,  resulting  in 
breast  deformation.  We  present  a  deformation  model  to  account  for  the  change  of  shape  and 
discuss  implications  in  image  reconstruction. 
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OCIS  codes:  (170.3830)  Mammography;  (170.6960)  Tomography;  (100.3010)  Image  reconstruction  techniques; 

(100.3190)  Inverse  problems 

1.  Introduction 

Near  Infrared  (NIR)  tomography  is  an  emerging  alternative  imaging  method  used  to  image  physiologic  parameters 
of  biological  tissue  in-vivo  such  as  water  and  hemoglobin  [1],  Measurements  of  light  propagation  (600  -  900  nm) 
within  tissue  can  be  used  to  map  internal  chromophore  concentrations  within  tissue.  Light  is  transmitted  through 
tissue  using  multiple  input  and  output  locations  on  the  surface  of  the  region  to  be  imaged  using  optical  fibers  for 
delivery  and  pickup  of  the  light  signals.  To  obtain  clinical  measurements  with  a  sufficient  signal  to  noise  ratio,  it  is 
key  to  ensure  good  contact  exists  between  the  optical  fibers  and  the  tissue.  The  breast  is  a  soft  tissue,  which  will 
deform  and  alter  its  shape  on  the  application  of  external  pressure.  The  amount  of  deformation  is  a  function  of  the 
tissue  mechanical  properties  and  the  amount  of  displacement  and  the  external  pressure  applied  by  the  optical  fiber 
array.  In  most  image  reconstruction  algorithms,  the  general  assumption  is  made  that  the  region  under  investigation  is 
a  uniform  circular  (2D)  or  conical  or  cylindrical  (3D)  domain  [2].  Little  work  has  been  done  to  evaluate  the  effect  of 
any  incorrect  geometry  in  image  reconstruction.  In  this  work,  the  effect  of  such  assumptions  is  investigated  by 
creating  a  deformation  of  the  breast  model  in  3D  to  generate  stimulated  clinical  data.  Images  are  then  reconstructed 
using  various  assumptions  regarding  the  imaging  domain  such  as  3D  conical  shaped  models  and  the  non-deformed 
breast  model  to  evaluate  image  quality. 

2.  Methods 

Soft  tissues  exhibit  nonlinear  elastic  behavior.  Nevertheless,  they  can  be  considered  a  linear  elastic  material  in 
situations  where  the  deforming  forces  produce  infinitesimal  deformations  (i.e.  <  5  %).  For  the  purpose  of  this  study, 
the  breast  was  modeled  as  linear  isotropic  pseudo-incompressible  medium  (i.e.  Poisson’s  ratio  (v)=  0.495).  Under 
these  assumptions  and  ignoring  internal  body  force,  the  governing  elasticity  equations  for  quasi-static  deformation  is 
given  by: 

( X  +  //)v(v  ♦  u)+  JuV2u  =  0  for  internal  nodes  in  domain  Q  (1) 

((/l  +  //)V(V  •  w)+  //V2w)l  h  =  h  for  nodes  on  the  boundary  8Q.  (2) 

Here  h  represents  a  unit  vector  directed  outwards  from  Q,  and  h  represents  the  traction  on  the  surface  or 
boundary  of  the  breast.  Note  that  u  represents  the  displacement  components  in  all  coordinate  directions,  and  \x  and 
X  are  Lame’s  elastic  constants.  For  an  isotropic  medium  these  constants  are  related  to  the  more  familiar  Young’s 

modulus  (E)  and  the  Poisson’s  ratio  (v)  by  — E —  and  X- - — - 

*  2(1  +  v)  (l  +  i/)(l  -  2v) 

It  should  be  stated  that  in  this  study  the  breast  was  assumed  to  be  traction  free  and  internal  body  forces  were 
neglected,  thus  the  problem  was  solved  by  imposing  a  prescribed  displacement  [3].  A  volume  mesh  of  a  female 
breast  of  a  volunteer  was  created  from  surface  image  data  that  was  acquired  using  a  3D  surface  camera  [Rainbow 
3D  Camera,  Genex  Technologies,  Kensington  MD].  The  3D  camera  projects  structural  illumination  patterns  onto 
the  object  and  calculates  3D  surface  profile  described  by  over  300,000  data  points  [4].  A  volume  mesh  was  then 
generated  using  the  Delaunay  algorithm.  The  mesh  has  a  geometry  of  130  x  136  x  60  mm,  Figure  1(a)  and  contained 
15501  nodes  corresponding  to  61 171  linear  tetrahedral  elements.  The  diameter  of  the  breast  at  its  mid-plane  where 
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optical  fiber  array  would  be  applied  is  approximately  88  mm.  To  calculate  the  deformation  due  to  16  equally  spaced 
optical  fibers  being  applied  at  the  mid-plane  of  the  breast,  i.e.  z  =  -30  mm,  it  is  assumed  that  each  optical  fiber 
pushed  the  breast  so  that  the  final  breast  diameter  at  z  =  -30  mm  is  70  mm,  and  that  the  diameter  of  each  optical  fiber 
is  6  mm.  The  modeled  elastic  properties  of  tissue,  were  assumed  as  isotropic  and  homogenous  with  Young’s 
Modulus  of  20  kPa  and  Poisson’s  ratio  of  0.495.  Further,  it  was  assumed  that  the  top  most  part  of  the  mesh,  i.e.  z  -  0 
mm  was  not  allowed  to  move  since  it  is  connected  to  the  chest.  Using  this  applied  displacement  as  a  boundary 
condition,  the  displacement  at  all  nodes  due  to  the  application  of  the  optical  fibers  was  calculated  and  a  deformed 
mesh  was  created,  Figure  1(b). 


Figure  1.  Volume  mesh  is  shown  of  the  (a)  normal  suspended  breast  and  (b)  the  deformed  mesh  after  the  application  of  the  optical  fiber  array. 

In  order  to  accurately  simulate  the  clinical  settings,  two  localized  anomaly  regions  were  placed  within  the  mid¬ 
plane  of  the  normal  breast,  both  at  z  =  -30  mm,  Figure  2(a).  First  was  an  absorbing  anomaly,  19  mm  from  the 
surface  and  a  radius  10  mm  with  a  coefficient  value  of  0.03  mm"1.  Second  was  a  reduced  scattering  anomaly,  also  19 
mm  from  the  surface  but  with  a  radius  of  7.5  mm  and  a  value  of  3  mm'1.  The  background  optical  properties  were 
modeled  with  an  absorption  coefficient  of  0.01  mm'1  and  a  reduced  scatter  coefficient  of  1.0  mm'1.  For  the  deformed 
breast  model  the  same  optical  properties  of  the  normal  breast,  including  the  anomalies  were  assumed.  Once  the 
deformation  of  the  model  was  calculated,  it  was  assumed  that  the  anomalies  were  free  to  move,  depending  on  the 
elastic  properties  of  the  breast,  Fjgure  2(b).  _  _____  _ 
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Figure  2.  2D  coronal  slices  through  (a)  the  normal  breast  mesh  and  (b)  the  deformed  breast  mesh,  showing  the  position  of  the  anomalies.  The 
right  hand  slice  is  near  the  chest  while  the  most  left  hand  slice  is  near  the  nipple. 

Using  the  deformed  mesh,  together  with  the  displaced  anomalies,  NIR  data  were  simulated  for  16  equally 
circularly  spaced  optical  fibers  placed  at  z  =  -30  mm.  Amplitude  and  phase  data  were  simulated  at  100  MHz,  and 
1%  noise  were  added.  This  data  was  then  used  as  simulated  patient  data  in  the  following  sections.  In  addition,  to 
allow  data  calibration  as  done  using  clinical  data,  and  discussed  elsewhere  [2],  the  data  was  simulated  for  16  equally 
spaced  optical  fibers  placed  in  a  circle  around  the  mid-plane  of  a  cylindrical  model  with  homogenous  background 
optical  properties. 

3.  Results 

Assuming  that  correct  anatomical  information  regarding  the  breast  is  available,  before  application  of  the  optical 
fiber  array,  the  normal,  non-deformed  mesh  was  used,  Figure  1(a),  during  reconstruction  to  produce  the  images 
shown  in  Figure  3(a).  Also,  images  reconstructed  assuming  a  conical  shaped  breast  (diameter  at  z  -  -30  set  to  equal 
the  diameter  of  the  optical  fiber  array)  are  shown  in  Figure  3(b).  Finally,  to  show  the  best  possible  reconstruction, 
the  deformed  mesh  was  used  for  reconstruction  of  images  and  these  results  are  shown  in  Figure  3(c). 

From  images  reconstructed  using  the  correct  but  non-deformed  breast,  although  the  absorbing  anomaly  has  been 
recovered,  the  reduced  scattering  images  contain  large  artifacts.  The  background  absorption  and  scatter  values  are 
about  0.006  mm'!  and  0.5  mm'1  respectively,  while  the  maximum  absorption  for  the  anomaly  is  0.014  mm"1.  In 
images  assuming  a  conical  model  the  background  absorption  and  scatter  values  are  about  0.007  mm'1  and  0.94  mm'1 
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respectively,  while  the  maximum  absorption  for  the  anomaly  is  0.014  mm'1.  Here,  although  the  target  values  are 
better  quantitatively,  the  images  contain  artifacts,  are  much  smaller,  and  the  absorption  scatter  cross-talk  is  more 
than  expected.  In  images  reconstructed  assuming  correct  3D  deformed  boundary,  both  the  absorbing  and  scattering 
regions  have  been  recovered  with  improved  localization  accuracy.  Some  cross  talk  exists  between  the  absorbing  and 
scattering  anomalies.  The  background  absorption  and  scatter  values  are  about  0.01  mm'1  and  0.98  mm'1  respectively, 
while  the  maximum  absorption  and  scatter  values  for  the  anomalies  are  0.013  mm’1  and  1.3  mm"1  respectively.  It  is 
important  to  note  that  although  the  images  in  Figure  3(a)  show  a  better  quantitative  accuracy  for  the  recovered 
absorbing  anomalies,  the  background  values  are  much  lower  than  expected  and  the  images  contain  far  more 
artifacts. 


Figure  3.  3D  simultaneous  reconstruction  of  absorption  and  reduced  scatter  from  deformed  breast  model  simulated  data  are  shown.  The  mesh 
used  for  the  reconstruction  was  (a)  the  normal  un-deformed  breast;  (b)  conical  mesh  and  (c)  the  correct  deformed  breast.  The  most  right  hand 
slice  is  near  the  chest  while  the  most  left  hand  slice  is  near  the  nipple  in  each  set  of  images. 


4.  Discussion 

In  the  case  where  it  is  assumed  that  the  exact  breast  geometry  is  known,  but  the  details  of  the  fibers  compression 
have  been  ignored,  although  the  absorption  anomaly  was  recovered  with  relatively  good  accuracy,  the  reduced 
scattering  image  contains  a  large  artifact.  Furthermore,  the  calculated  background  values  for  absorption  and  scatter 
are  0.006mm'1  and  0.5  mm'1  respectively,  which  are  much  lower  than  expected.  These  images  are,  perhaps,  not  as 
accurate  and  useful  as  the  conical  geometry  case,  since,  although  there  is  a  more  accurate  3D  model,  the  diameter  of 
the  breast  was  not  correct,  particularly  within  the  measurement  plane,  i.e.  z  =  -30  mm.  When  using  the  geometrically 
correct  deformed  mesh  to  reconstruct  images  of  absorption  and  scatter,  superior  localization  accuracy  resulted.  The 
scattering  object  has  also  been  recovered  in  this  case.  However,  the  quantitative  accuracy  of  the  recovered  anomalies 
is  not  as  good  with  maximum  values  for  absorption  and  scatter  of  0.013mm"1  and  1.3  mm’1,  respectively.  This 
modest  degradation  in  quantitative  accuracy  has  been  reported  elsewhere  and  is  a  common  problem  in  3D  imaging 
algorithms,  which  may  be  solved  through  the  use  of  multi-stage  algorithms  [5]  or  with  inclusion  of  a  priori  data. 
Nonetheless  using  the  correct  model  with  correct  information  regarding  the  deformed  boundary  has  produced 
images  with  very  little  to  no  artifacts. 
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Abstract 

A  multi-spectral  direct  chomophore  and  scattering  reconstruction  technique  has  been 
implemented  for  near-infrared  frequency-domain  tomography  in  recovering  images  of  total 
hemoglobin,  oxygen  saturation,  water  and  scatter  parameters.  The  method  applies  the  spectral- 
constraint  of  the  chromophores  and  scattering  spectra  into  the  inversion  algorithm,  thereby 
reducing  the  parameter-space  of  the  inversion  process.  This  new  method  was  validated  in 
experimental  phantoms  and  clinical  images.  Analysis  of  the  recovered  images  indicates  that  the 
added  spectral-constraint  in  the  inversion  allows  a  lower  level  of  regularization  to  be  used  to 
achieve  the  same  minimization  error  norm  in  the  objective  function.  This  approach  also  produces 
a  significant  reduction  in  image  artifacts  and  is  improved  in  image  accuracy. 


1 


Near  Infra-Red  (NIR)  tomography  is  known  for  its  ability  to  distinguish  between 
malignant  and  normal  tissue  based  on  a  high-contrast  mechanism  arising  from  intrinsic  processes 
such  as  angiogenesis  and  hypoxia.  Absorption-based  parameters  can  be  recovered  such  as  total 
hemoglobin,  oxygen  saturation  and  water  fraction  as  well  as  elastic  scattering.  These  can  provide 
measures  of  physiology  or  pathophysiology  related  to  the  intravascular  and  extravascular  spaces 
as  well  as  the  composition  of  tissue.  In-vivo  studies  have  demonstrated  levels  of  hemoglobin  in 
tumors  over  twice  that  in  normal  breast  [1,  2];  however,  one  of  the  current  challenges  is  to 
optimize  the  quantitative  accuracy  with  which  these  parameters  can  be  determined.  Typically, 
spectral  tomography  approaches  make  use  of  sparse  data  at  discrete  wavelengths  instead  of  a 
complete  spectrum.  This  sparse  spectral  sampling  coupled  with  an  image  reconstruction  process 
which  is  ill-posed,  amplifies  errors  in  quantifying  the  physiological  parameters  of  the  tissue. 

The  method  illustrated  here  reduces  noise  in  the  recovery  of  chromophores,  especially 
significant  in  water  and  scattering,  by  eliminating  the  need  to  reconstruct  the  absorption  and 
scattering  coefficients  as  intermediate  endpoints.  Instead,  the  chromophore  concentrations  and 
scatter  parameters  are  estimated  directly  by  incorporating  the  known  Beer’s  law  relation  and  Mie 
scattering  behavior  as  constraints.  Corfu  et  al  [3]  have  implemented  this  approach  using 
continuous  wave  measurements  to  find  the  optimal  four  wavelengths  that  reduce  the  crosstalk 
between  absorption  and  scatter  parameters.  Their  results  from  simulations  are  encouraging  under 
the  assumption  that  all  change  in  scattering  is  due  to  scatter  amplitude  with  scatter  power  kept 
constant.  A  similar  approach  to  CW  data  has  been  implemented  by  Li  et  al[4]  who  have  used  two 
of  three  wavelengths  under  the  assumption  that  there  is  no  scattering  perturbation.  They  have 
applied  this  method  to  find  chromophore  concentrations  directly,  and  have  shown  in  simulated 
and  experimental  data  that  their  technique  results  in  reduced  artifacts  and  crosstalk.  In  this  study, 
we  extend  the  approach  to  application  of  full  frequency-domain  data  sets  at  six  wavelengths,  and 
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demonstrate  reconstruction  in  phantoms  and  humans.  The  method  uses  a  finite  element  model  of 
the  diffusion  equation  and  reconstructs  images  for  five  parameters:  oxyhemoglobin,  hemoglobin, 
water  fraction,  scatter  amplitude  and  scatter  power,  with  no  assumptions  on  scatter  amplitude  or 
power.  The  results,  in  addition  to  showing  a  reduction  in  noise  in  the  recovered  chromophore 
concentrations,  also  provide  evidence  of  improvement  in  image  accuracy  and  suppression  of 
image  artifacts. 

The  NIR  frequency  domain  imaging  system  (tailored  for  breast  imaging)  has  been 
described  in  detail  elsewhere  [5].  Briefly,  it  consists  of  optical  fibers  placed  in  three  planes  in  a 
circular  geometry.  Each  plane  has  16  source-detector  positions  and  intensity  modulated  light  at 
100MHz  is  used  at  six  different  wavelengths  in  the  range  660-850nm.  The  signal  amplitude  and 
phase  is  calibrated  to  compensate  for  system  offsets  by  matching  measured  data  from 
homogeneous  calibration  phantoms  to  simulated  results  from  the  finite  element  model. 

Image  reconstruction  uses  Newton’s  method  to  iteratively  solve  for  the  optical  properties 
where  the  underlying  model  is  a  finite  element  solution  of  the  diffusion  approximation  to  the 
radiative  transfer  equation  and  appropriate  boundary  conditions  have  been  applied  (type  Ill-Robin 
type  boundary).  Further  details  can  be  found  in  [6],  In  simple  terms,  this  iterative  procedure 
involves  the  Taylor  series  approximation: 

<3^  =  35//  ^ 

where  ^  refers  to  the  change  in  boundary  data,  ^  is  the  Jacobian,  the  matrix  containing  the 
sensitivity  of  the  boundary  data  to  change  in  optical  properties,  3  =  [3/(  ,3^] ;  and  is  the 
update  in  the  optical  properties  given  by  dju  =  (d/ua,dfc)  where  f.ia  is  the  absorption  coefficient 
in  mm'1  and  K  is  the  diffusion  coefficient  given  by  K  = - \ - in  mm,  and  jus  is  the  reduced 

3(//„+/0 

scattering  coefficient  in  mm'1.  The  inverse  problem  corresponding  to  equation  (1)  is  solved  with 
application  of  the  Levenberg-Marquardt  regularization  scheme  for  stabilization.  Previously,  the 
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optical  properties  at  each  wavelength  were  obtained,  and  then  the  calculation  of  the  chromophore 
concentrations  was  performed  in  a  single  step  using  the  linear  relation 

(2) 


where  e  is  the  molar  absorption  spectra  of  the  tissue  absorbing  chromophores,  in  our  case,  oxy¬ 
hemoglobin  (Hb02),  hemoglobin  (Hb)  and  water;  and  c  is  the  concentration  of  these 
chromophores.  Similarly,  for  scattering  the  approximation  to  Mie  theory  [7,  8], 


was  used  to  derive  scatter  amplitude  (a)  and  scatter  power  (b)  with  wavelength  in  pm.  Assuming 

~  _j¥  „  _  f¥ 

that  we  know  and  as  calculated  by  the  previous  method,  in  the  new  approach 

the  measurements  at  all  wavelengths  are  coupled  together,  and  the  relationships  in  equations  2 
and  3  are  combined  to  create  a  new  set  of  equations,  which  for  each  wavelength,  is  represented 
by: 


d<f>x  —  *sc  xdc  +  Ada  +  -3b  Adb 


(4) 


where 


^c.K  — 


8<j>  |  d</>  dju 


dc  d/J,  dc  ^  for  each  chromophore  in  the  model.  From  equation  2,  we  get 

dpi 

df*  ~  £^c ,  so  that,  substituting  for  dc  > 


81 

dc 


_d$_  ,  [d£ 

d/u  x  dju 


(5) 


where  <8>  refers  to  the  Kronecker  tensor  product. 

Similarly,  expressions  for  A  and  A  can  be  derived,  and  come  out  to  be 


o  ,  = 
a, A 


da  8k  da 


(6) 
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Similarly,  for  the  scatter  power  b  part: 

3m  =  3r  (-  3k2  )(-  In  A)\  A  (7) 

The  overall  system  of  equations  is  built  by  substituting  the  relations  from  (5),  (6)  and  (7)  into 
equation  4: 


_ 

fdc  ) 

(t>K ) 

^  c** 

^  C\M"*  c2M^  c3,/ll^  aj\^  Ml 

dc2 

— 

C**  C**  C** 
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•• 
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The  size  of  the  left  hand  vector  is  number  of  wavelengths  times  number  of  measurements  per 
wavelength  while  the  size  of  the  right  hand  side  vector  is  number  of  chromophores  times  number 
of  nodes  for  updating  each  chromophore  in  the  mesh.  The  individual  Jacobians  for  each 
chromophore  were  computed  using  a  dual  mesh  technique  [9],  on  a  2000  node  mesh  for  forward 
diffusion  calculations  and  a  425  node  mesh  for  image  reconstruction.  Equation  8  was 
implemented  by  building  up  the  new  Jacobian  (5x425  by  480x6).  The  same  Levenberg- 
Marquardt  regularization  scheme  was  applied  and  the  computation  time  was  approximately 
25minutes  for  typically  5-7  iterations;  the  measure  of  convergence  being  when  the  difference 
between  measured  and  calculated  data  (projection  error)  was  less  than  2%  of  previous  iteration 
value. 

In  order  to  test  the  accuracy  of  the  reconstruction  procedure  with  experimental 
measurements,  data  was  collected  on  a  liquid  tissue-simulating  phantom  within  a  plastic  circular 
container  consisting  of  1%  Intralipid  and  1%  blood  in  buffered  saline.  The  blood  hematocrit  was 
measured  before  the  experiment  by  a  clinical  co-oximeter  that  showed  1%  whole  blood  was 
equivalent  to  17pM  of  hemoglobin;  and  water  is  expected  to  be  100%.  Measurements  were 
carried  out  on  this  phantom  in  fully  oxygenated  and  deoxygenated  states,  with  the  latter  done  by 
reducing  pC>2  in  the  phantom  through  the  addition  of  yeast.  Figure  1  shows  the  mean  and  standard 
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deviation  from  recovered  images  for  total  hemoglobin  (Hbr),  oxygen  saturation  (St02),  water 
fraction  and  scatter  parameters  by  both  spectrally-constrained  approach  and  the  former  method 
where  each  wavelength  is  reconstructed  separately.  The  mean  recovered  HbT,  S,02  and  water 
fraction  values  are  accurate  with  less  than  3%  mean  error  for  spectral  method  and  5.1%  for  the 
separate  wavelength  reconstruction  using  the  stopping  criterion  of  less  than  2%  projection  error 
change  between  iterations.  The  error  bars  in  the  figure  show  the  standard  deviation  in  the 
homogeneous  images;  the  standard  deviations  in  the  spectrally-constrained  images  are  much 
lower  than  those  from  the  conventional  technique.  The  images  were  not  spatially-filtered  in  either 
method,  so  that  no  averaging  of  properties  occurred  in  order  to  preserve  maximum  quantitative 
accuracy.  The  reconstruction  by  spectral  method  converges  at  iteration  6,  and  is  -12  iterations  for 
conventional  technique,  for  each  wavelength.  Reduced  standard  deviation  in  the  old  technique 
can  be  achieved  by  artificially  stopping  the  reconstruction  before  the  convergence  criterion;  the 
standard  deviation  using  iteration  6  for  conventional  technique  is  shown  in  Figure  1,  and  is  lower 
than  obtained  by  letting  the  conventional  reconstruction  converge  till  projection  error  changes  by 
less  than  2%.  Reduced  standard  deviation  in  the  conventional  technique  can  also  be  achieved  by 
the  use  of  filtering,  with  the  tradeoff  in  quantitative  damping  of  optical  properties;  whereas  the 
spectral  technique  offers  lower  noise  in  images  without  using  filters.  The  reduction  in  standard 
deviation  is  especially  significant  in  water  fraction  and  scattering  parameters  and  artifacts 
occurring  close  to  the  boundary  in  the  images  have  been  suppressed  using  the  spectral  approach. 
The  spectral  reconstruction  technique  was  also  applied  to  the  complete  set  of  16  measurement 
sets  obtained  for  varying  p02  measured  separately  by  a  chemical  microelectrode.  The  recovered 
oxygen  saturation  follows  the  predicted  Hill  curve  with  a  mean  error  of  7.7%  (data  not  shown). 

To  evaluate  the  spectrally-constrained  method  with  heterogeneous  experimental  data,  a 
similar  liquid  phantom  as  above  was  constructed  with  1%  intralipid  and  1%  pig  blood  (measured 
to  be  -0.01 7mM)  in  buffered  saline  in  a  90mm  diameter  plastic  bottle.  A  clear  plastic  cylindrical 
inclusion  of  diameter  20mm  was  placed  at  10  o’clock  in  the  phantom  10mm  from  the  edge  with 
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the  same  Intralipid  solution  but  twice  the  concentration  of  blood  (2%  =  0.034  mM).  Amplitude 
and  phase  measurements  were  obtained  in  the  frequency  domain  at  six  wavelengths  and  the 
images  were  reconstructed  by  applying  both  the  spectral  and  conventional  two-stage  method,  as 
shown  in  Figure  2.  Recovered  images  from  the  direct  approach  are  shown  in  the  top  row  and  were 
obtained  with  the  stopping  criterion  of  less  than  2%  error  change  between  iterations  (7  iterations). 
This  technique  converges  faster  than  the  conventional  two-stage  method.  The  corresponding 
images  obtained  by  using  absorption  and  scatter  estimates  at  the  seventh  iteration  from  the 
conventional  method  are  shown  in  middle  row  of  Figure  2  while  those  images  based  on  the  same 
projection  error  change  criterion  for  each  wavelength  appear  in  the  bottom  row.  The  direct 
chromophore  reconstruction  shows  quantitatively  accurate  property  recoveiy  with  8%  error  in 
HbT  and  less  than  2%  error  in  St02  and  water  fraction,  in  the  region  of  interest.  The  images  in  the 
top  and  middle  row  are  comparable;  however,  the  images  from  the  middle  row  show  an  artifact 
developing  in  the  center,  which  is  more  pronounced  in  the  bottom  row  when  the  same  stopping 
criterion  is  applied.  The  images  are  expected  to  be  homogeneous  in  all  other  parameters  except 
total  hemoglobin.  The  water  fraction  and  scattering  images  in  particular  show  lower  noise  in  the 
spectral  reconstruction.  Some  inhomogeneity  is  observed  in  the  water  fraction  images,  which 
results  from  the  walls  of  the  plastic  inclusion.  No  filtering  has  been  applied  to  any  of  the  images, 
and  some  artifacts  are  observable  close  to  the  edges  in  certain  images,  which  may  be  due  to  the 
boundary  of  the  plastic  phantom.  In  general,  the  coupling  of  all  measured  wavelengths  to 
reconstruct  chromophores  and  scattering  directly  appears  to  suppress  the  common  noise 
components  within  the  individual  datasets,  leading  to  more  artifact  free  images. 

The  direct  spectral  reconstruction  method  has  also  been  applied  to  clinical  tomography 
data  obtained  from  measurements  on  a  66  year  old  subject  with  a  5-10  mm  spiculated  mass 
evident  on  mammography,  which  was  later  diagnosed  as  Invasive  Ductal  Carcinoma.  The  subject 
underwent  the  NIR  exam  in  accordance  with  the  Dartmouth  protocol,  where  informed  consent 
was  obtained  prior  to  the  exam.  The  position  of  the  anomaly  was  provided  through 
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mammography  and  determined  to  be  at  1 1:30  o’clock  in  the  cranocaudal  view  and  located  about 
lcm  from  surface.  NIR  data  was  collected  in  three  planes.  Reconstructed  images  are  shown  in 
Figure  3  for  the  plane  in  line  with  the  tumor.  A  localized  increase  in  HbT  was  observed  and  the 
contrast  available  was  1.7: 1.0  in  tumor  versus  background.  The  St02  image  showed  a  decrease  at 
the  location  of  the  tumor  with  a  contrast  of  0.7:1.  The  corresponding  St02  image  with  the 
conventional  technique  (not  shown)  was  noisier,  and  a  similar  decrease  was  not  observed.  The 
water  fraction  image  showed  heterogeneity  in  the  range  0.25  to  0.5%,  which  is  considerably 
smaller  than  the  range  of  0  to  1%  observed  with  the  old  method.  Tighter  data  ranges  are  similarly 
observed  in  the  scatter  images,  and  artifacts  evident  in  the  separate  wavelength  technique  have 
been  suppressed  in  the  spectral  approach. 

In  summary,  the  new  direct  spectral  reconstruction  technique  has  been  applied  to 
homogeneous  and  heterogeneous  data  and  the  results  have  consistently  reduced  artifacts  as  well 
as  improved  quantitative  accuracy  in  recovered  parameters.  With  this  new  technique,  higher 
hemoglobin  content  at  the  site  of  tumor  is  observed  in  the  clinical  case,  as  well  as  lower  oxygen 
saturation.  The  water  fraction  and  scattering  images  show  the  most  significant  improvement  in  all 
cases  by  reduction  of  the  image  artifacts. 


This  work  has  been  supported  through  NIH  grants  PO1CA80139  and  R01CA69544  and 
Department  of  Defense  DAMD1 7-03-01-0405.  The  authors  would  like  to  acknowledge  expert 
clinical  collaboration  with  Steven  P.Poplack,  Sandra  K.  Soho  and  Christine  Kogel. 
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Figure  Captions 

Figure  1.  Mean  and  standard  deviations  from  images  of  HbT  [mM],  S,C>2  [%],  water  fraction, 
lOOxa  and  lOOxb  for  the  spectral  approach  using  convergence  criterion  of  projection  error 
variation  less  than  2%  between  iterations  (converges  at  iteration  6);  for  the  conventional  method 
using  corresponding  6,h  iteration  values  and  for  the  conventional  method  using  the  convergence 
criterion  as  applied  for  spectral  method  (converges  at  ~12  iterations  for  each  wavelength).  These 
were  generated  from  experimental  measurements  on  a  homogeneous  phantom  containing  17pM 
whole  blood  with  1%  Intralipid  in  saline. 

Figure  2.  Reconstructed  images  of  HbT  [mM],  S,02  [%],  water  fraction,  ‘a’  and  ‘b’  parameters 
from  experimental  data  on  a  phantom  solution  of  Intralipid  in  saline  with  0.017  mM  hemoglobin 
in  background  and  0.034  mM  in  the  heterogeneity.  Images  with  the  spectral  approach  are  shown 
in  row  1  using  a  stopping  criterion  of  projection  error  variation  being  less  than  2%  (~7  iterations), 
Row  2  shows  the  conventional  method  after  7  iterations,  and  row  3  shows  the  old  method  using 
the  same  stopping  criterion  as  row  1  (10-12  iterations  for  each  wavelength). 

Figure  3.  HbT  [mM],  S,02  [%],  water  fraction,  ‘a’  and  ‘b’  images  are  shown  from  measurements 
on  a  cancer  subject  with  a  5-10  mm  Invasive  Ductal  Carcinoma  at  the  1 1 :30  clock  face  position. 
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ABSTRACT 


The  design  and  implementation  of  a  multi-spectral,  frequency-domain  near 
infrared  tomography  system  is  outlined,  which  operates  in  an  MRI  magnet  for 
utilization  of  MR-guided  image  reconstruction  of  tissue  optical  properties.  Using 
long  silica  optical  fiber  bundles,  measurements  of  light  transmission  through  up  to 
12  cm  of  female  breast  tissue  can  be  acquired  simultaneously  with  MRI  scans. 
The  NIR  system  utilizes  six  optical  wavelengths  from  660-850  nm  using  intensity 
modulated  diode  lasers  nominally  working  at  100  MHz.  Photomultiplier  tube 
detector  gain  levels  are  electronically  controlled  on  a  time  scale  of  200  ms, 
thereby  allowing  rapid  switching  of  the  source  to  locations  around  the  tissue. 
There  are  no  moving  parts  in  the  detection  channels  and  for  each  source  position, 
15  PMTs  operating  in  parallel  allow  sensitivity  down  to  0.5  pW/cm2  at  the  tissue 
surface.  Images  of  breast  tissue  optical  absorption  and  reduced  scattering 
coefficients  are  obtained  using  a  Newton-type  reconstruction  algorithm  to  solve 
for  an  optimal  solution  using  the  measurement  data.  In  medical  imaging,  it  is 
beneficial  to  compare  the  same  tissue  volume  as  seen  by  a  variety  of  modalities, 
and  perhaps  more  importantly,  there  is  the  hypothesis  that  one  imaging  system 
which  has  high  spatial  resolution  can  be  used  to  enhance  the  reconstruction  of 
another  system  which  has  good  contrast  resolution.  In  this  study  we  explore  the 
synergistic  benefits  of  a  combined  NIR-MRI  data  set,  specifically  the  ways  in 
which  MRI  (i.e.  high  spatial  resolution)  enhances  NIR  (i.e.  high  contrast 
resolution)  image  reconstruction.  The  design,  calibration,  and  performance  of  the 
imaging  system  are  described  in  the  context  of  preliminary  phantom  tests  and 
initial  in  vivo  patient  imaging. 


L  INTRODUCTION 


Diffuse  optical  tomography  (DOT)  with  near  infrared  (NIR)  light  can  be  used  to  produce 
spatially  resolved  images  of  tissue  optical  properties.  These  optical  property  maps  can  be 
acquired  at  different  wavelengths  and  combined  to  reveal  hemoglobin  concentration,  oxygen 
saturation,  water  and  fat  content,  as  well  as  a  description  of  scattering  structures1-5.  These  latter 
parameters  are  important  indicators  of  metabolic  activity,  functional  processes,  or  presence  and 
staging  of  disease.  NIR  diffuse  tomography  generally  suffers  from  comparatively  low  spatial 
resolution  due  to  the  multiple  scattering  events  that  occur  along  each  photon  path6.  However,  the 
promise  of  this  imaging  modality  lies  in  the  fact  that  it  affords  new  physical  bases  for  contrast  in 
tissue.  For  example,  hemoglobin-based  contrast  in  tumors  relative  to  normal  tissue  is 
exceptionally  high  (i.e.  between  100%-300%).  However,  improving  the  limitations  associated 
with  its  low  spatial  resolution  is  fundamental  to  implementing  this  technology  clinically.  A 
priori  knowledge  of  tissue  structure  can  be  used  to  constrain/guide  the  iterative  NIR  image 
reconstruction  process,  and  improve  the  spatial  resolution  and  quantitative  accuracy  of  recovered 
physiological  parameters.  Consequently,  NIR  techniques  have  been  combined  with  several  high 
spatial-resolution,  structure-bearing  imaging  modalities  including  x-ray  tomosynthesis7, 
ultrasound  (US)8,  and  magnetic  resonance  imaging  (MRI)9'11,  to  study  human  tissues  and  small 
animals.  In  this  report,  we  present  a  combined  NIR-MRI  system  for  imaging  female  breast 
tissue,  and  explore  the  benefits  of  the  combined  data  set  in  a  planar  tomographic  geometry  where 
the  breast  is  imaged  pendent  in  a  standard  MR  breast  coil. 

In  recent  years,  several  important  technological  factors  have  contributed  to  the 
advancement  of  NIR  spectral  imaging  applied  to  breast  cancer  characterization,  including  an 
increased  understanding  of  light  and  tissue  interaction,  and  new  software-based  developments  in 
image-reconstruction  algorithms  for  DOT5.  NIR  radiation  (700-900  nm)  is  non-ionizing,  light 
delivery  and  detection  instrumentation  is  relatively  inexpensive,  and  uncomfortable  tissue 
compression  is  not  required  (as  in  mammographic  level  compression).  NIR  light  transmission 
can  typically  be  measured  through  as  many  as  12  cm  of  breast  tissue  depending  on  its 
composition  and  the  wavelength  of  the  NIR  source.  Photon  propagation  within  the  breast  is  well 
described  by  diffusion  theory  since  the  probability  of  photon  scattering  is  much  greater  than 
absorption.  Light  transmission  measurements  can  be  combined  with  diffusion  theory  to  provide 
robust  image  reconstruction  of  tissue  optical  coefficients.  We  use  a  Newton-type  algorithm  to 
solve  for  the  optimal  solution  that  provides  a  minimum  error  between  the  measured  data  and 
predicted  response  from  a  model  of  the  frequency-domain  diffusion  equation12.  Frequency- 
domain  measurements  contain  information  on  both  the  amplitude  change  in  the  optical  flux  and 
time-based  information.  In  theory,  the  approach  allows  separation  of  the  absorption  and  reduced 
scattering  coefficients  (pa  and  p/)  because  the  signal  amplitude  and  phase  shift  provide  a  unique 
data  set  to  solve  the  estimation  problem  associated  with  recovering  both  coefficients 
simultaneously. 

A  variety  of  imaging  methods  have  achieved  high-spatial-resolution  imaging  with 
acceptable  to  excellent  soft  tissue  contrast,  including  X-ray  computed  tomography  (CT),  US,  and 
MRI.  These  techniques  primarily  provide  images  of  tissue  structure  and  have  a  limited  ability  to 
monitor  parameters  related  to  tissue  function  other  than  through  the  introduction  of  exogeneous 
contrast  agents.  Alternatively,  nuclear  medicine  approaches  are  routinely  used  to  image  tissue 
functions,  such  as  metabolic  fluorodeoxyglucose  uptake13,  and  many  commercial  systems  exist 
to  co-register  these  images  with  the  structural  data  derived  from  CT,  US  and  MRI.  The 


combination  of  high  resolution  structural  imaging  with  lower  resolution  functional  information  is 
a  major  emphasis  in  contemporary  medical  imaging,  and  customized  hybrid  imaging  systems  are 
being  developed  to  avoid  the  complications  associated  with  tissue  movement  between  imaging 
exams,  which  compromises  the  accuracy  of  post-procedure  co-registration. 

The  application  of  NIR  tomography  to  provide  spatially-resolved  functional  information, 
such  as  hemoglobin  levels,  oxygen  saturation,  water,  lipid  and  scatterer  content  will  likely  be 
important,  yet  customized  imaging  systems  which  couple  to  MRI,  US  or  CT  need  to  be 
developed  to  evaluate  and  exploit  this  potential.  The  pioneering  work  of  Ntziachristos  et  al.](> 
was  the  first  to  demonstrate  feasibility  of  the  hybrid  NIR-MR  approach  for  clinical  breast 
imaging.  Their  approach  consisted  of  a  co-planar  array  of  optical  source/sensors  detecting  time- 
resolved  illuminations  when  applied  to  the  surface  of  the  breast.  Data  was  presented  as  maps  of 
localized  optical  response  at  each  detector  site  but  was  otherwise  not  processed  into  depth- 
resolved  images  of  optical  property  distributions.  Zhu  et  al ,8  have  demonstrated  clinical 
application  of  a  combined  NIR-US  system.  With  a  hand-held  probe,  they  simultaneously 
measure  reflection  of  ultrasound  and  modulated  NIR  excitations  at  the  breast  tissue  surface,  and 
construct  co-registered  images  of  structure  and  optical  absorption  and  hemoglobin  concentration. 

One  of  the  most  challenging  issues  in  hybrid  system  development  is  the  incorporation  of 
anatomical  information  as  prior  constraints  into  NIR  image  reconstruction.  This  has  been 
explored  theoretically  in  several  successful  reports14'19.  Structural  information  can  be  used  to 
guide  functional  image  reconstruction  through  knowledge  of  tissue  composition  and  location  by: 
(i)  alteration  of  the  objective  functions  used  in  image  reconstruction7’20,  (ii)  reduction  of  the 
number  of  unknown  parameters  by  treating  regions  of  the  same  tissue  type  as  single  zones10,21, 
or  (iii)  introduction  of  special  regularization  schemes  that  can  stabilize  the  inverse  problem  and 
emphasize  image  contrasts7,22’23.  Optimizing  spatial  resolution  more  likely  will  depend  on  the 
application  of  all  three  techniques,  while  also  being  strongly  influenced  by  the  signal  to  noise 
ratio  of  the  measurements,  the  optical  contrast  available,  and  the  number  of  projections  used. 
Some  investigators  have  successfully  combined  different  approaches  with  multi-stage  image 
reconstruction  algorithms17’ 18.  However,  there  is  still  no  clear  consensus  on  how  best  to  utilize 
the  structural  information  to  enhance  or  improve  the  recovery  of  functional  NIR  information. 

This  work  describes  the  first  system  to  combine  multi-spectral  frequency-domain  NIR  in 
a  planar  tomographic  geometry  with  MRI  for  imaging  breast  tissue.  NIR  tomography  has  shown 
the  ability  to  localize  changes  in  functional  tissue  parameters  in  vivo,  and  MRI  has  the  advantage 
of  offering  a  particularly  rich  amount  of  anatomical  information,  specifically  about  the  layered 
adipose  and  glandular  tissue  structure  of  the  breast.  This  system  is  designed  to  combine  the 
benefits  of  both  modalities  into  the  construction  of  a  single  quantitative  image.  In  particular, 
since  NIR  tomography  has  significant  difficulty  with  the  recovery  of  properties  in  layered 
structures,  the  initial  information  from  MRI  can  significantly  improve  the  estimation  of  breast 
properties,  including  in  localized  regions  such  as  tumors. 

The  design  and  operation  of  the  NIR-MRI  system  elements  are  described  in  Section  II. 
The  NIR  component  is  similar  to  an  imaging  system  described  previously24,  which  is  currently 
being  evaluated  clinically,  yet  has  the  unique  design  feature  of  having  no  moving  parts  in  the 
detection  channels,  allowing  significant  reduction  in  the  NIR  data  acquisition  time.  In  Section 
III  we  outline  the  theoretical  basis,  and  the  practical  application  and  utilization  of  image 
reconstruction  in  NIR  tomography.  After  applying  several  approaches  to  optimizing  this  hybrid 
reconstruction,  an  algorithm  is  examined  that  takes  advantage  of  the  composite  data  set.  In 


Section  IV  we  discuss  system  performance  and  present  images  of  phantom  and  breast  optical 
properties  which  are  both  high  resolution  and  quantitatively  accurate. 


II.  SYSTEM  DESIGN 

This  section  describes  the  four  elements  that  comprise  the  NIR-MRI  system:  (A)  light 
delivery,  (B)  detector  array,  (C)  fiberoptic  patient  interface,  and  (D)  computer  control  and 
electronics.  Element  (C)  extends  from  (A)  and  (B)  into  an  MRI  scanner  for  clinical  studies,  as 
shown  schematically  in  Figure  1.  Photographs  of  the  rack-mounted  system  and  patient  interface 
are  presented  in  Figure  2.  Section  (E)  briefly  describes  phantom  fabrication,  and  its  importance 
in  imaging  system  development  and  validation. 

A.  Light  delivery 

The  system  deploys  six  laser  diodes:  661  nm  (40  mW),  752  nm  (50  mW),  785  nm  (50  mW),  805 
nm  (50  mW),  829  nm  (50  mW),  and  849  nm  (50  mW).  Each  wavelength  is  amplitude  modulated 
at  100  MHz  by  mixing  a  DC  current  source  (LDX-3220,  ILX  Lightwave,  Bozeman,  MT)  and  an 
AC  current  from  a  frequency  generator  (IFR-2023A,  IFR  Systems,  Wilmington,  MA),  through  a 
bias  T  (#5545,  Picosecond  Pulse  Labs,  Boulder,  CO).  Each  diode  is  held  in  a  laser  tube 
(Thorlabs,  Newton,  NJ),  and  mounted  on  a  linear  translation  stage  (MA2515P5-S2.5  Velmex, 
Bloomfield,  NY).  This  stage  directs  a  specified  wavelength  into  one  of  sixteen  bifurcated  optical 
fiber  bundles  which  were  custom  designed  for  this  application  (Ceramoptec,  East  Longmeadow, 
MA).  The  248-piece  bundles  (0.37  N.A.,  0.68  packing  fraction)  are  pure  silica  core  (210  pm), 
silicone  clad  (230  pm)  fibers  suitable  for  transmission  wavelengths  from  400  nm  to  2400  nm. 

The  source  light  is  delivered  through  the  central  seven  fibers  in  each  bundle,  and  the  remaining 
fibers  surrounding  these  are  delivered  to  the  detectors.  The  common  end,  which  makes  contact 
with  the  tissue,  has  a  diameter  of  4mm.  Each  fiber  bundle  is  13  meters  in  length  and  extends 
from  the  instrument  cart,  located  outside  of  the  MR  suite,  into  the  bore  of  the  scanner  (1.5T 
whole  body  imager,  GE  Medical  Systems,  Milwaukee,  WI)  to  the  patient  interface.  The 
efficiency  of  the  optical  switching  is  approximately  50%,  yielding  an  average  source  power  of  15 
mW  at  the  tissue  surface. 

B.  Light  detection 

For  each  source  excitation,  light  transmission  is  recorded  from  15  surface  locations.  This  signal 
is  measured  by  15  photomultiplier  tubes  (PMT  R6357,  Hamamatsu,  Japan)  operating  in  parallel. 
The  gain  of  the  PMTs  is  varied  to  account  for  the  large  variation  in  light  level  between  detectors 
depending  on  their  distance  from  the  source.  The  gains  are  set  by  applying  computer  generated 
voltages  between  0.4  and  1.2  V  to  their  control  lines,  which  sets  the  anode  to  cathode  voltage 
between  approximately  350  V  to  1000  V,  respectively.  Using  the  higher  gain  settings,  a  PMT 
can  reliably  measure  optical  signals  in  the  pW  range.  The  optimal  gain  levels  are  determined 
prior  to  each  imaging  series.  Each  PMT  is  fixed  to  a  particular  fiber,  so  it  is  necessary  to  switch 
gains  electronically  during  the  course  of  data  collection.  A  100  MD  resistor  was  used  in  the 
dynode  chain  of  each  PMT  to  achieve  fast  settling  times  after  gain  adjustment  (200  ms  for  large 
gain  changes).  Electrical  heterodyning  through  RF  mixers  (Minicircuits,  Brooklyn,  NY)  is  used 
to  down  convert  the  100  MHz  PMT  signal  to  a  lower  frequency  (500  kHz).  This  offset 
frequency  is  achieved  with  a  second  frequency-synthesizer  which  is  synchronized  to  the  one 


driving  the  laser  current,  and  is  set  to  100.0005  MHz.  The  resulting  offset  frequency  is  filtered 
and  amplified  by  a  16  channel  circuit  designed  for  this  application  (Audon  Electronics, 
Nottingham,  UK),  then  read  by  the  computer.  Lock-in  detection  is  executed  in  software  to 
extract  amplitude  and  phase  data  for  each  of  the  detectors  in  parallel. 

C.  Fiberoptic  patient  interface 

The  MR  exam  is  performed  using  a  breast  array  coil  (MRI  Devices,  Waukesha,  WI)  that  offers 
high-resolution  imaging.  The  coil  also  provides  an  open  architecture,  which  allows  for  the 
integration  of  the  NIR-breast  interface  shown  in  Figure  1.  A  circular  ring  machined  from 
polyvinyl  chloride  (PVC)  positions  the  common  ends  of  the  sixteen  fibers  around  the  full 
circumference  of  the  pendant  breast.  Phosphor  bronze  compression  springs  (k=0.091bs/in,  Ace 
Wire  Spring  &  Form,  McKees  Rocks,  PA)  guide  each  fiber  through  holes  in  the  ring,  into  light 
contact  with  the  tissue  surface.  The  ring  separates  into  equal  halves  so  that  it  can  easily  be 
moved  from  one  breast  to  the  other.  The  angular  separation  between  each  fiber  is  21  degrees, 
except  between  fibers  adjacent  to  the  line  of  ring  devision,  where  the  separation  is  33  degrees. 
The  ring  can  be  positioned  vertically,  such  that  the  plane  of  measurements  intersects  the  region 
of  interest  in  tissue.  This  design  allows  each  fiber  to  move  independently,  therefore  their  radial 
positions  may  vary  by  several  millimeters.  In  order  to  avoid  serious  artifacts  in  reconstructed 
images,  the  position  of  each  fiber  must  be  accurately  known  during  data  acquisition.  Annular 
fiducial  markers  (MM  3005,  IZI  Medical  Products,  Baltimore,  MD)  are  fixed  to  each  fiber  and 
can  be  located  with  sub-millimeter  accuracy  in  the  MRI. 

D.  Computer  system 

A  PC  running  Labview  software  (National  Instruments,  Austin,  TX)  is  used  to  control  all  light 
delivery  and  detection  equipment.  The  laser  current  source  and  frequency  generator  parameters 
are  set  by  a  general  purpose  interface  bus  (GPIB,  NI).  The  linear  translation  stage  is  addressed 
through  the  serial  port.  An  analog  output  board  (NI)  is  used  for  PMT  gain  control.  A 
multipurpose  data  acquisition  (DAQ)  board  acquires  the  16  analog  input  channels  and  the  single 
reference  channel.  This  board  also  provides  six  digital  output  lines  to  the  high  power  radio¬ 
frequency  switch  for  the  laser  sources.  For  each  source  position,  15  signals  from  the  detector 
system  are  amplified  by  a  gain  of  1000  and  low  pass  filtered  to  prevent  aliasing  prior  to  the  DAQ 
board  using  a  16  channel  amplifier  and  filter  network  mounted  in  a  BNC  coupled  box  (Audon 
Electronics,  Nottingham  UK).  Data  are  acquired  for  500  ms,  and  phase  and  amplitude  of  each 
signal  are  calculated  and  written  to  file.  Including  the  time  required  to  determine  optimal  gain 
values,  measurement  with  6  wavelengths  takes  approximately  4  min.  The  MR  exam  is 
controlled  separately,  operated  in  parallel,  and  a  full  volume  breast  MRI  is  of  similar  duration.  A 
FORTRAN,  or  MATLAB  based  reconstruction  program  reads  and  processes  the  NIR  data25'27. 
MR  images  are  processed  off-line  with  an  addition  to  the  MATLAB  software  package,  and 
incorporated  into  a  modified  iterative  optical  property  reconstruction  (see  Section  III). 

E.  Phantom  Design 

Tissue-simulating  phantoms  with  known  property  distributions,  geometries,  and  imaging 
orientations  are  commonly  used  to  validate  imaging  systems.  We  have  developed  a  recipe  for 
producing  gelatin  materials  with  desired  optical  properties,  and  a  shelf  life  of  several  months28. 

A  heated  mixture  of  water,  gelatin  (G2625,  Sigma  Inc.),  India  ink  (for  absorption),  and  titanium 
dioxide  powder  (for  scatter)  (Ti02,  Sigma  Inc.)  is  poured  into  a  mold  of  a  desired  shape,  and 


solidified  by  cooling  to  room  temperature.  Variation  in  the  water  concentration  provides  MR 
contrast,  and  variable  gel  stiffness.  The  phantom  imaged  here  (Section  IV.  B.)  combines  three 
gels  with  different  optical  properties  in  an  irregular  structure. 

HI.  DATA  PROCESSING  AND  IMAGE  RECONSTRUCTION 

Quantitative  NIR  imaging  with  model-based  methods  requires  (A)  important  instrument 
calibration  procedures,  and  (B)  a  reconstruction  algorithm  that  incorporates  an  accurate  model  of 
light  propagation  in  tissue. 

A.  System  calibration 

Calibration  issues  and  other  practical  considerations  associated  with  our  NIR  imaging  approach 
have  been  discussed  in  detail  elsewhere26.  Two  important  procedures  are  briefly  noted  here:  (1) 
detector  calibration,  and  (2)  homogeneous  phantom  calibration.  First,  the  amplitude  and  phase 
response  of  each  detection  channel  must  be  characterized  in  order  to  remove  systematic  noise  in 
the  data  acquisition  hardware.  Each  detector  is  exposed  to  the  same  optical  signal,  and  the 
differences  in  log  amplitude  and  phase  are  used  as  correction  factors.  The  log  amplitude 
response  of  the  PMT  is  plotted  against  the  log  of  the  input  power  for  each  gain  setting.  A  log — 
log  regression  is  performed  and  the  coefficients  are  used  to  calibrate  detected  PMT  amplitude  in 
terms  of  optical  power.  The  phase  does  not  fluctuate  significantly  with  changing  light  level  for  a 
single  gain  setting  (i.e.  minimal  phase-amplitude  cross-talk),  but  is  altered  dramatically  with 
changing  gain.  Relative  phase  differences  between  detectors  are  stored  for  calibration.  These 
calibration  curves  are  very  similar  to  those  created  by  McBride  et  a/24.  This  characterization 
needs  to  be  performed  only  once  as  long  as  the  system  is  not  modified. 

The  second  important  practical  procedure  is  the  correction  for  inter  fiber  variations  and 
coupling  issues,  which  is  accomplished  through  a  homogeneous  phantom  calibration  process26,29. 
This  accounts  for  offsets  due  to  optical  fiber  differences  in  transmission  and  alignment,  as  well 
as  any  errors  in  discretization  or  data-model  mismatch.  A  homogeneous  phantom  is  generally 
measured  each  day,  and  after  system  changes.  The  differences  between  data  measured  from  the 
phantom,  and  data  calculated  from  the  model  are  stored  and  subtracted  from  measurements  of 
the  heterogeneous  phantom  or  tissue  under  investigation.  A  homogeneous  fitting  algorithm  is 
used  to  determine  the  pa  and  p/  values  supplied  to  the  model  calculation.  This  algorithm  can 
also  be  used  to  calculate  the  initial  optical  properties  specified  in  iterative  reconstruction  of 
heterogeneous  media.  When  dealing  with  tissues  having  arbitrary  shape,  the  effectiveness  of  this 
fitting  algorithm  and  homogeneous  phantom  calibration  hinges  on  the  accurate  specification  of 
source  and  detector  locations.  The  ability  to  extract  accurate  fiber  positions  from  MRI  scans 
preserves  the  integrity  of  this  method  for  non-uniform  boundary  data. 

B.  FEM  analysis 

Data  acquired  from  the  detection  system  is  processed  by  a  FEM  based  reconstruction  algorithm 
to  generate  tomographic  images  of  absorption  and  reduced  scattering  coefficients  simultaneously. 
The  algorithm  exploits  the  frequency-domain  diffusion  equation  approximation  to  light  behavior 
in  a  highly  scattering  medium25, 

-V  •  Z)(r)VO(r,<9)  +  //»  +  —  Wr.ffl)  =  S(r,co)  (1) 

v  c  J 

where  S(r,co )  is  an  isotropic  light  source  at  position  r,  @(r,co)  is  photon  density  at  r,  c  is  the 


speed  of  light  in  tissue,  a>  is  the  frequency  of  light  modulation,  //„  is  the  absorption  coefficient, 
and  D  =  - - - — —  is  the  diffusion  coefficient.  The  reduced  (transport)  scattering  coefficient  is 

X “a+Ms) 

given  by  ju's  =  fis  (1  -  g) ,  where  g  is  the  mean  cosine  of  the  single  scatter  function  (the 
anisotropy  factor),  and  fis  is  the  scattering  coefficient . 

For  a  given  pa  and  p/  distribution,  the  diffusion  equation  is  used  to  predict  the  optical 
flux  at  the  detector  sites  for  each  source  excitation.  In  the  inverse  problem  (image 
reconstruction),  the  goal  is  the  estimation  of  optical  properties  at  each  FEM  node,  based  on 
measurements  of  optical  flux  at  the  detector  sites  on  the  tissue  surface.  This  is  achieved 
numerically  by  minimizing  the  difference  between  the  calculated  data  <PC ,  and  measured  data, 
0M ,  for  all  source/detector  combinations  (NM).  Typically, 

NM  .  x, 

=!(<-<?  (2) 

1=1 

is  minimized  in  a  least  squares  sense  by  setting  the  first  derivative  equal  to  zero,  and  using  a 
Newton-Raphson  approach  to  find  the  set  of  optical  property  values  which  approximate  the  point 
of  stationarity.  We  use  a  Levenberg  Marquardt  algorithm,  and  repeatedly  solve  the  equation 

a  =  (jTJ  +  AlYjTb,  (3) 

where  b  =  (<£>r  -  &M )?  is  our  data  vector,  and  a  is  the  solution  update  vector,  a  =  [SD] ;  SjuaJ  j, 

defining  the  difference  between  the  true  and  estimated  optical  properties  at  each  reconstructed 
node  j.  Here,  A  is  a  regularization  factor  to  stabilize  matrix  inversion  and  J  is  the  Jacobian 
matrix  for  our  model,  which  is  calculated  using  the  Adjoint  method30. 

Improving  NIR  reconstructions  by  incorporating  MRI  data  have  been  explored  in 
previous  work19  ,  and  by  other  authors10, 18,21.  Techniques  used  to  produce  the  images  shown 
in  Section  IV  include:  (i)  accurately  defining  the  imaging  volume,  (ii)  tailoring  a  regularization 
scheme  which  optimizes  the  reconstructed  contrast  of  a  suspicious  area  in  the  image,  and  (iii) 
reducing  the  number  of  unknown  parameters  by  segmenting  tissue  types  visible  in  the  MRI. 
Defining  the  imaging  volume,  (i),  is  relatively  straight  forward  and  can  be  accomplished  by 
creating  a  structured  finite  element  mesh  from  the  MRI.  This  mesh  will  often  contain  regional 
differences  depending  on  tissue  types  present,  and  an  irregular  outer  boundary  due  to 
impressions  caused  by  fiber-tissue  contact.  In  simulation  studies,  not  presented  here,  it  has  been 
found  that  image  reconstruction  accuracy  is  easily  degraded  if  the  mesh  (2D  or  3D)  does  not 
represent  the  true  outer  tissue  boundary.  Step  (ii)  can  be  accomplished  by  changing  Al  in  Eq.  3  to 
XA,  where  A  is  a  regularization  matrix,  or  filter  matrix.  Regularization  can  be  thought  of  as  a 
smoothing  operator,  where  one  can  apply  selective  smoothing  by  linking  together  the  property 
updates  for  all  nodes  associated  with  the  same  region  or  tissue  type.  Modifications  to  the 
Jacobian  matrix  size  in  a  parameter  reduction  technique  are  used  to  implement  step  (iii)20. 


IV.  PERFORMANCE  AND  EXPERIMENTAL  RESULTS 

In  this  section  we  (A)  compare  measurements  from  this  system  to  those  from  a  similar 
NIR  tomography  instrument  described  previously24  for  a  number  of  tissue  simulating  phantoms. 
NIR-MRI  phantom  studies  are  described  in  (B),  and  in  vivo  images  are  shown  in  (C).  For  all 


results  presented  here,  we  used  a  single  laser  diode  (785  nm)  for  convenience,  and  2D  modeling 
and  image  reconstruction.  System  performance  and  image  quality  at  the  other  five  wavelengths 
are  comparable,  and  3D  imaging  is  readily  achievable. 

A.  System  performance 

Measurement  repeatability  in  terms  of  phase  and  log  amplitude  error  was  assessed  by  serially 
imaging  a  phantom  with  optical  properties  similar  to  those  of  average  breast  tissue  (pa=0.004, 
Ps/=1 .35).  The  average  RMS  error  at  each  detector  site  was  determined  to  be  0.26%  in  ac 
intensity  and  1 .04°  in  phase.  These  values  compare  with  those  obtained  from  the  system 
described  by  McBride  et  al.  (0.32%  in  ac  intensity  and  0.48°  in  phase)24.  RMS  error  for  each  of 
the  240  source-detector  combinations  for  both  systems  is  plotted  versus  PMT  signal  in  Figure  3. 
For  both  phase  and  amplitude,  error  sharply  increases  when  incident  light  falls  below 
approximately  0.5  pW.  These  points  are  excluded  in  the  calculation  of  average  RMS  error.  As 
an  additional  comparison,  we  used  both  devices  to  measure  a  collection  of  phantoms  (N=15)  of 
varying  diameters  and  optical  properties.  After  processing  the  measurements  with  the 
homogeneous  phantom  calibration  procedure  described  in  Section  III,  we  used  the  homogeneous 
fitting  algorithm  to  determine  a  single  pa  and  p/  for  each  material.  Table  1  shows  the  optical 
coefficients  obtained  with  the  NIR-MRI  system,  along  with  a  measure  of  their  discrepancy  with 
those  obtained  with  our  stand-alone  DOT  instrumentation.  We  observed  good  agreement 
between  the  two  systems.  The  absorption  and  scattering  coefficients  show  correlation 
coefficients  of  R2=0.984  and  R2=0.980,  respectively. 

B.  Phantom  Imaging 

Due  to  the  limited  spatial  resolution  of  DOT,  layered  media,  small  objects,  and  low  contrast 
heterogeneity  pose  key  challenges  in  image  reconstruction.  The  capability  of  the  presented 
system  to  address  these  challenges  was  investigated  by  imaging  a  phantom  comprised  of  three 
gels  with  different  optical  properties.  The  phantom  is  cylindrical  and  the  boundary  between  the 
outer  layer  and  inner  layer  is  irregular.  A  two  centimeter  spherical  inclusion  is  embedded  within 
the  inner  layer.  The  optical  coefficients  of  each  gel  are  known  (Table  2),  as  each  material  was 
created  using  a  practiced  recipe  to  give  desired  values.  Furthermore,  the  true  value  of  each  layer 
was  validated  by  creating  a  separate  homogenous  phantom  for  each  layer  (at  the  same  time  as 
creating  layered  phantom)  and  measuring  bulk  properties  with  a  homogeneous  fitting  algorithm. 
To  increase  MRI  contrast  between  layers,  Omniscan™  (gadodiamide)  was  added  to  the  inner 
layer  (0.005  g/ml). 

A  photograph  of  the  phantom,  alongside  another  spherical  inclusion  is  shown  in  Figure 
4(a).  Figure  4(b)  shows  a  T1 -weighted,  gradient  echo  MRI  (25  ms  TR,  3  ms  TE,  45°  flip  angle) 
cross-section  of  the  phantom  in  the  plane  of  the  optical  fibers.  This  was  used  to  create  a  2D 
structured  finite  element  mesh  and  to  locate  the  positions  of  the  16  fibers  (Figure  4(c)).  Fiber 
fiducials  were  not  used  in  this  experiment,  but  the  16  impressions  caused  by  each  optical  fiber 
are  clearly  visible  around  the  perimeter  of  the  gel.  Each  of  the  three  regions  are  also  visible, 
corresponding  to  each  of  the  three  types  of  optically  variant  gel.  3D  meshing  and  imaging  is  also 
possible,  given  that  a  stack  of  MR  slices  represent  the  full  volume  structure  of  the  phantom.  The 
system  design  allows  for  NIR  data  acquisition  in  multiple  planes  if  desired. 

Using  log  amplitude  and  phase  data,  images  of  optical  properties  were  reconstructed  with 
two  different  algorithms.  The  first  algorithm  solves  Eq.  (3)  without  a  priori  guidance,  except  for 
the  use  of  the  mesh  and  source  locations  from  Figure  4(c).  The  corresponding  reconstructed 


images  are  presented  in  Figure  5(a).  The  second  algorithm  uses  the  layered  structure  from  the 
MRI  to  constrain  Eq.  (3).  The  regularization  parameter  associated  with  each  reconstructed  node 
is  adjusted  based  upon  its  location  and  area  of  influence.  Lower  values  are  used  for  small 
regions  close  to  the  center  of  the  reconstructed  volume,  whereas  peripheral  regions  (which  are 
prone  to  large  artifacts)  have  larger  regularization.  A  regularization  matrix  (A),  or  filter  matrix, 
based  upon  the  a  priori  structural  information  from  the  MRI  was  also  used  to  further  improve  the 
algorithm.  The  optimal  values  for  X  and  A  were  determined  in  simulation  studies  of  this 
geometry,  using  similar  contrast  and  noise  levels.  These  can  also  be  chosen  automatically  once 
an  empirical  knowledge  of  their  effect  is  established.  The  images  reconstructed  from  this 
modified  constrained  algorithm  are  shown  in  Figure  5(b). 

Qualitatively,  the  algorithm  that  uses  MRI  information  to  guide  the  iterative  process 
performs  much  better.  The  optical  property  images  in  Figure  5(a)  are  blurry  and  “edge  artifacts” 
are  clearly  visible,  especially  for  p/.  The  absorbing  sphere  near  the  center  of  the  phantom  is 
recovered  with  poor  spatial  resolution,  and  its  contrast  is  underestimated.  The  property  maps  in 
Figure  5(b)  have  improved  resolution.  Spatial  resolution  and  contrast  of  the  spherical  inclusion 
are  better.  As  a  quantitative  measure  of  the  accuracy  of  image  reconstruction,  we  compute  the 
RMS  error  between  the  target  image  and  the  reconstructions,  node  by  node.  This  measure 
indicates  a  dramatic  improvement  for  the  second  method  [RMS(pa)=0.0024,  RMS(p/)=0.19] 
relative  to  the  first  [RMS(pa)-0.0034,  RMS(ps/)=0.25]. 

C.  In  Vivo  Imaging 

The  Institutional  Review  Board  (IRB)  at  the  Dartmouth  Hitchcock  Medical  Center  approved  this 
clinical  examination  protocol,  and  written  informed  consent  was  obtained  from  all  volunteering 
women.  Figure  6  (a  and  b)  show  a  photograph  of  a  subject  lying  on  the  examination  platform, 
and  an  anatomical  axial  MR  image  through  the  breast.  NIR  measurements  (4  wavelengths  in  this 
case)  and  a  full  volume  MRI  (50  coronal  slices,  25  ms  TR,  6  ms  TE,  45°  flip  angle,  2  mm  slice 
thickness)  were  obtained  in  less  than  10  minutes.  Figure  7(a)  shows  an  anatomically  coronal 
MRI  with  16  fiducial-marked  fibers  (appearing  as  bright  white  spots  outside  the  breast).  The 
radiographic  density  of  this  participant  is  heterogeneously  dense  (HD),  and  a  large  interior  region 
of  glandular  tissue  is  easily  defined  in  the  FEM  mesh  shown  in  Figure  7(b). 

As  with  the  phantom  study,  we  present  images  of  pa  and  p/  at  785nm  reconstructed  with 
two  different  algorithms.  The  first  result,  obtained  by  solving  Eq.  (3),  without  a  priori  guidance, 
is  shown  in  Figure  7(c).  As  expected,  we  see  higher  absorption  and  scatter  in  the  glandular 
region  (central)  relative  to  the  adipose  tissue  (peripheral).  However,  the  region  of  increase  does 
not  span  the  full  area  expected,  and  heterogeneity  is  visible  (especially  around  the  perimeter  of 
the  image.  The  second  algorithm  assumes  homogeneous  optical  properties  for  each  tissue  type, 
and  utilizes  parameter  reduction20,  which  leads  to  a  ‘fitting’  for  four  values:  Pa  adipose-0.003,  ps7 
adipose=0.93,  Pagianduiar=0’-006,  ps/g!anduiar=l-12  (Figure  7(d)).  This  algorithm  is  robust  to  noise  and 
converges  after  a  few  iterations.  The  result  is  quantitatively  logical,  and  similar  values  are 
recovered  using  spatially  varying  regularization  (as  described  in  the  phantom  reconstruction). 


V.  DISCUSSION  AND  CONCLUSION 


This  article  describes  a  novel  in  vivo  breast  imaging  system  that  synergistically  combines 
NIR  tomography  with  MRI.  Tissue  structures  visible  with  high  resolution  in  MRI  can  be  applied 


a  priori  to  optical  property  reconstructions  from  frequency-domain  NIR  measurements.  Thus, 
the  reconstruction  process  can  be  optimized  to  produce  high  resolution,  quantitatively  accurate 
maps  of  absorption  and  reduced  scattering  coefficients,  and  ultimately  physiologically  relevant 
parameters.  Various  physiological  tissue  types  exhibit  significant  contrast  in  the  NIR,  primarily 
though,  the  combined  system  could  be  very  effective  at  locating  and  diagnosing  breast  tumors. 

The  NIR  component  provides  multi-spectral  (6  wavelengths)  frequency-domain  (log 
amplitude  and  phase)  data  from  16  fiber-tissue  contact  positions  around  the  breast’s 
circumference.  Data  acquisition  is  fully  automated,  and  a  complete  set  of  measurements  (240 
source-detector  pairs)  for  one  wavelength  requires  approximately  40  s.  Light  detection  is  highly 
sensitive  (subpicowatt  limit)  with  low  noise  (RMS  <0.26%,  <1.04°  in  phase).  All  optical 
elements  and  controls  are  mounted  in  a  portable  cart,  and  operation  inside  strong  magnetic  fields 
is  facilitated  with  long  optical  fibers  and  an  easy-to-use  positioning  system/patient  interface. 

The  deployment  of  dual  modality  NIR  imaging  systems  in  clinical  applications  has  been 
limited  to  date,  mainly  because  of  the  complexity  of  the  image  reconstruction  problem.  Here,  a 
representative  phantom  and  in  vivo  study  using  one  wavelength  (785  nm)  are  presented.  We 
have  shown  that  co-registered  MRI  validates  and  improves  optical  property  estimation  in  2D 
tomographic  image  reconstructions  when  specialized  algorithms  are  used.  Future  work  will 
involve  3D  modeling  and  reconstruction,  which  could  further  improve  both  qualitative  and 
quantitative  aspects  of  the  recovered  coefficient  values. 

Preliminary  results  are  encouraging,  and  have  allowed  us  to  optimize  reconstruction 
techniques,  and  automate  constraint  selection.  We  have  performed  several  phantom  studies,  and 
demonstrated  the  feasibility  of  imaging  volunteers  with  healthy  breasts.  Through  the  study  of 
more  healthy  subjects,  with  different  radiographic  densities,  we  aim  to  compare  functional 
parameters  of  adipose  versus  glandular  tissue.  To  test  the  system’s  ability  at  diagnosing  tumors, 
cancer  patients  will  be  recruited,  and  the  simultaneous  exam  will  likely  involve  MR  contrast 
enhancement. 
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6  LDs  on  translation  stage  A16  bifurcated  optical  fibers  (13  m) 


Figure  1.  Schematic  design  of  a  dual  modality  NIR-MRI  system  (left).  Frequency-domain  NIR 
tomography  is  performed  inside  the  MRI  unit  (upper  right).  Six  laser  diodes  (660-850nm)  are 
amplitude  modulated,  and  sixteen  projections  yield  240  measurements  of  amplitude  and  phase  of 
transmitted  light.  The  fiber  optic  array  is  positioned  inside  the  open  breast  coil,  to  allow 
positioning  along  the  length  of  the  pendant  breast  (bottom  right). 


Figure  2.  Photograph  of  the  rack  mounted,  portable  system  (left).  System  components  are 
marked,  including  (A)  frequency  generators,  (B)  optical  switching  stage,  (C)  PMT  detection 
plate,  and  (D)  optical  fibers  and  patient  interface.  The  fiber-patient  interface  (at  right)  can 


accommodate  breasts  6-12  cm  in  diameter.  A  close  up  of  the  tissue  coupling  system  shows  that 
the  fibers  are  spring-loaded  and  make  light  contact  with  the  full  circumference  of  a  pendant 
breast. 


Figure  3.  Repeatability  assessment  for  (a)  log  amplitude  and  (b)  phase  of  NIR  component  of 
combined  NIR-MRI  system  and  stand-alone  NIR  tomography  system  at  Dartmouth.  The 
performance  of  the  two  systems  is  comparable.  In  routine  operation,  PMT  voltage  signals  are 
above  1.0x1  O'5  V. 


Ha  (mm'1) 
(NIR-MRI) 

%  difference 

p/  (mm1) 
(NIR-MRI) 

%  difference 

Mean 

0.0055 

-6.2 

1.34 

9.1 

Max 

0.0102 

6.3 

1.91 

29.6 

Std.  Dev. 

0.0024 

8.7 

0.41 

7.9 

Table  1.  Optical  properties  at  785nm  for  15  tissue  simulating  phantoms,  and  %  differences  from 
values  determined  using  another  NIR  tomography  system. 


Outer  Layer 

Inner  Layer 

Inclusion 

Ha  (mm1) 

H»;  (mm1) 

0.0044 

0.0062 

0.02 

0.51 

0.68 

0.9 

Table  2.  Approximate  optical  properties  of  gelatin  phantom. 


Figure  4.  (a)  The  gelatin  phantom,  (b)  T1 -weighted  MRI  image  of  the  structure  of  the  phantom 
cross-section.  The  inner  layer  (light  region)  contains  gadodiamide  for  MR  contrast  (0.005  g/ml). 
The  inner  inclusion  is  a  2  cm  diameter  sphere,  (c)  Finite  element  mesh  derived  from  the  three¬ 
layered  structure  image  in  (b).  The  circles  on  the  outer  boundary  indicate  the  fiber  locations. 
Axes  are  in  millimeters. 


Figure  5.  Images  in  (a)  show  optical  property  reconstructions  without  MRI  guidance.  In  (b), 
reconstructions  improve  when  the  interior  structural  information  of  the  MRI  is  incorporated. 
Comparing  these  two  images,  the  latter  has  a  reduction  in  artifacts  in  the  reduced  scattering 
image,  and  both  the  spatial  resolution  and  contrast  have  improved.  The  recovered  values 
compare  favorably  with  the  approximate  true  values  shown  in  Table  2. 


Figure  6.  Photograph  in  (a)  shows  a  female  volunteer  prepared  for  the  simultaneous  MRI-NIR 
exam.  In  an  anatomically  axial  T1  -weighted  MRI  slice  from  the  right  breast,  (b),  fiducial 
markers  (outside  the  breast)  indicate  the  location  of  the  fiber  plane. 


(c)  *10  (d)  *10 

Figure  7.  In  (a),  an  anatomically  coronal  T1 -weighted  MRI  displays  adipose  (outer)  and 
glandular  (inner)  tissue  types.  A  two  layer  structured  FEM  mesh  and  source  locations,  created 
from  (a),  is  shown  in  (b).  Optical  property  reconstructions  in  (c)  are  obtained  without  utilizing 
the  internal  structure  of  (b).  In  (d),  MRI  data  guides  a  two-region  parameter  fitting  algorithm. 
This  result  is  improved  relative  to  (c),  and  shows  increased  absorption  and  scatter  in  glandular 
relative  to  adipose  tissue. 


